{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as stats\n",
    "\n",
    "import pandas as pd \n",
    "import seaborn as sns\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Correlated Time Series\n",
    "\n",
    "## Review\n",
    "\n",
    "We have been modeling stock prices of stock $k$ as Geometric Brownian Motion, \n",
    "\n",
    "$$\\frac{dS_k}{S_k} = \\mu dt + \\sigma dW_k,$$\n",
    "\n",
    "where $\\langle dW_k\\rangle = 0, \\langle dW_k^2\\rangle = dt$.\n",
    "We then included correlation, meaning that the error terms of two stocks may not be independent of each other, as originally assumed:\n",
    "\n",
    "$$\\langle dW_1 dW_2\\rangle = \\rho_{1,2} dt\\ne 0.$$\n",
    "\n",
    "$\\rho_{1,2}$ is called the correlation coefficient betwen 1 and 2.\n",
    "\n",
    "We proceeded to define the covariance matrix which we could extend to $N$ stocks and allow for more generality, that is, $\\sigma_i\\ne \\sigma_j$ for $i\\ne j$.\n",
    "\n",
    "$$\\Sigma = \n",
    "\\begin{pmatrix}\n",
    "\\sigma_1^2 & \\rho_{1,2} \\sigma_1\\sigma_2 & \\dots & \\rho_{1,N} \\sigma_1\\sigma_N  \\\\ \n",
    "\\dots & & & \\dots \\\\\n",
    "\\rho_{N,1}\\sigma_N\\sigma_1  & \\dots &\\dots & \\sigma_N^2\n",
    "\\end{pmatrix}.$$\n",
    "\n",
    "It is convenient to write this in terms of the covariance, $\\sigma_{1,2} = \\rho_{1,2}\\sigma_1\\sigma_2$.\n",
    "\n",
    "## Single Factor Model\n",
    "\n",
    "In today's lecture, we have shifted our focus a little bit:\n",
    "We are considering a time series for excess returns that can be described as\n",
    "\n",
    "$$R_{i,t}^e = \\beta_i F_t + \\sigma_i \\epsilon_{i,t},$$\n",
    "\n",
    "where $F_t$ is a time series of the \"returns\" of one risk factor that all stocks are exposed to in some way. We write $\\langle F_t\\rangle =\\lambda$ and $\\langle F_t^2\\rangle = \\sigma_F^2$.\n",
    "Their sensitivity to the risk factor is given by their factor loading $\\beta_i$. \n",
    "We still have the random error term, and by removing the one common risk factor, we assume $\\epsilon_{i,t}$ is now truly idiosyncratic.\n",
    "\n",
    "What is the result of this? \n",
    "Now the covariance matrix is \n",
    "\n",
    "$$\\Sigma = \n",
    "\\begin{pmatrix}\n",
    "b_1^2 \\sigma_F^2 +\\sigma_1^2 & b_1 b_2 \\sigma_F^2 & \\dots & b_1 b_N \\sigma_F^2  \\\\ \n",
    "\\dots & & & \\dots \\\\\n",
    "b_N b_1\\sigma_F^2  & \\dots &\\dots & b_N^2 \\sigma_F^2 + \\sigma_N^2\n",
    "\\end{pmatrix}.$$\n",
    "\n",
    "#### Problem 1.\n",
    "Let's run a little experiment. Let's build three groups of stocks, $N_1=10, N_2=N_3=5$ where the stocks in the first group have $b=0.8$, in the second group $b=0.3$, and in the third group $b=0$.\n",
    "\n",
    "Assume that $\\lambda = 8\\times 10^{-4}$, $\\sigma_F=3\\times 10^{-2}$, and $\\sigma_i=2\\times10^{-2}$ for all $i$. Generate 100 data points and plot the correlation matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD4CAYAAAAZ1BptAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOy9eZQkV30u+N3Ycs/aq3pfpN7U2lFLgA2SgGchDot4BsbSYJDH2LIZM9jjGT/jsYFnbL9jv+czfocz2B4ZsMEzbCPbSBjZGCyBwQhJLam1tNS7Wt3V3dW1V+6ZEZF3/oi4ETciIyIjqzKzqrrud06frsqKzIqsjLi/+33fbyGUUggICAgICARBWu0TEBAQEBBYuxBBQkBAQEAgFCJICAgICAiEQgQJAQEBAYFQiCAhICAgIBAKZbVPoJsYHR2lu3btWu3TEBAQEFhXeOaZZ2YppWNBP7uigsSuXbtw+PDh1T4NAQEBgXUFQshrYT8TcpOAgICAQCi6EiQIIXcTQo4TQk4RQj4R8PPbCSHPEkIMQsj7fT+7nxBy0v53P/f4LYSQF+3X/CwhhHTjXAUEBAQE4mPFQYIQIgP4HIB3ADgI4D5CyEHfYecA/AKAr/ieOwzg0wBeD+A2AJ8mhAzZP/4LAA8A2Gv/u3ul5yogICAg0Bm6wSRuA3CKUnqGUtoA8DUA9/AHUErPUkpfAND0PfftAL5LKZ2nlC4A+C6AuwkhmwHkKaVPUKtvyJcBvLcL5yogICAg0AG6ESS2AjjPfT9pP7aS5261v277moSQBwghhwkhh2dmZmKftICAgIBAe3QjSAR5BXG7BoY9N/ZrUkofpJQeopQeGhsLzOASEBAQEFgmuhEkJgFs577fBuDiCp87aX+9nNcUEBAQEOgSuhEkngawlxCymxCiAbgXwCMxn/sdAHcRQoZsw/ouAN+hlF4CUCSEvMHOavowgIe7cK4CAgICDnSziW8cPo9mU4xMCMOKgwSl1ADwMVgL/isAvkEpPUoI+Qwh5D0AQAi5lRAyCeADAP5vQshR+7nzAP4AVqB5GsBn7McA4KMAPg/gFIDTAP5ppecqICAgwOPHp+fwnx56AS9cWFrtU1mz6ErFNaX0UQCP+h77FPf10/DKR/xxXwTwxYDHDwO4rhvnJyAgIBCEmm4CAKoNc5XPZO1CVFwLCAhsWOimlZXfMP3Z+QIMIkgICAhsWBim5UXUdcEkwiCChICAwIZFQzCJthBBQkBAYMOCMYmGIYJEGESQEBAQ2LBgnkRdBIlQiCAhICCwYeEY1yJIhEIECQEBgQ0LXchNbSGChICAwIaFKzeJ7KYwiCAhICCwYSHkpvYQQUJAQGDDgslNdZECGwoRJAQEBDYsBJNoDxEkBAQENixECmx7iCAhICCwYSGym9pDBAkBAYENCyE3tYcIEgICAhsWIgW2PUSQEBAQ2LAQvZvaQwQJAQGBDQvRBbY9RJAQEBDYsHDkJl0EiTCIICEgILBh4chNgkmEoitBghByNyHkOCHkFCHkEwE/TxBCvm7//ElCyC778Q8SQo5w/5qEkJvsn33ffk32s/FunKuAgIAAQ0NkN7XFioMEIUQG8DkA7wBwEMB9hJCDvsM+AmCBUroHwJ8B+BMAoJT+v5TSmyilNwH4EICzlNIj3PM+yH5OKZ1e6bkKCAgI8DBEkGiLbjCJ2wCcopSeoZQ2AHwNwD2+Y+4B8CX764cAvI0QQnzH3Afgq104HwEBAYFYcHo3iSARim4Eia0AznPfT9qPBR5DKTUALAEY8R3zc2gNEn9tS02fDAgqAgICAiuCaMvRHt0IEkGLN+3kGELI6wFUKKUvcT//IKX0egBvtv99KPCXE/IAIeQwIeTwzMxMZ2cuICCwoeFWXItiujB0I0hMAtjOfb8NwMWwYwghCoABAPPcz++Fj0VQSi/Y/xcBfAWWrNUCSumDlNJDlNJDY2NjK3gbAgICGw1CbmqPbgSJpwHsJYTsJoRosBb8R3zHPALgfvvr9wN4jFJKAYAQIgH4ACwvA/ZjCiFk1P5aBfAuAC9BQEBAoIvQuWI6e0kS8EFZ6QtQSg1CyMcAfAeADOCLlNKjhJDPADhMKX0EwBcA/C0h5BQsBnEv9xK3A5iklJ7hHksA+I4dIGQA3wPwVys9VwEBAQEejElQChhNClUW1qcfKw4SAEApfRTAo77HPsV9XYPFFoKe+30Ab/A9VgZwSzfOTUBAQCAMOldEVzeaUGVRX+yH+IsICAhsWOhmE5odGEStRDBEkFgFLFV1lOrGap+GgMCGh2FSpBMyABEkwiCCxCrgf/nqc/g//v7F1T4NAYENDUopGmYT2YSluosgEYyueBICnWGmWEeppq/2aQgIbGgYTcu0ZkFCDB4KhggSq4CGYTo9YwQEBFYHrANsWrPkJlErEQwRJFYBuklRN4QnISCwmmAdYLNJ1fO9gBfCk1gF6GYThaoIEgICqwmW/pq1jWsxeCgYIkisAnSziapuenK0BQQE+gtXbrKNa3E/BkIEiVUAy6Io1gSbEBBYLbhMQmQ3RUEEiVUAawVQqIoMJwGB1YI/SIjspmCIILEKYBenYBICAqsHtlkTxXTREEGiz2g2qZOfXdhAtRLn5yu478GfYKZYX+1TERAAIOSmuBBBos/Qm+6FuJHkppcuLOGJM3P41vP+USMCAqsDFiQywriOhAgSfQajuMDGYhKsUOmfX5pa5TMRELDA7sUM8yRECmwgRJDoM3SO0m4kT4KZgk+/Ni8kJ4E1gRa5STCJQIgg0WfwtREbSW5iTIJS4F9eFmxCYPXB7sWUaMsRCREk+gx+t1LYSEzCpvLjuYSQnATWBJjclFAkaIokUmBDIIJEn7FxPQnrBnzPjVvwxOk5LFYaq3xGAhsdjEkoMkFClkR2UwhEkOgzvHLTBmISRhOEAO+6cQuMJsX3Xple7VMSWAf4xb95Gn/x/dM9eW12L6qyhIQqgkQYRJDoM/gLsVdMYqHcwDs/+0OcmSn15PWXg7rRREKRcOO2AWwZSArJSSAWnj23gBcvLPbktRmr12QJmiwJTyIEXQkShJC7CSHHCSGnCCGfCPh5ghDydfvnTxJCdtmP7yKEVAkhR+x/f8k95xZCyIv2cz5LCCHdONfVhrt7IT3LbjozW8LRiwW8dLHQk9dfDuq6iYQigxCCt1+3Cf92ckaMcBWIBKUUxZrRs/uEl5s0RTCJMKw4SBBCZACfA/AOAAcB3EcIOeg77CMAFiilewD8GYA/4X52mlJ6k/3vV7nH/wLAAwD22v/uXum5rgWw3ctwRutZdlOlYen/5TW0CDMmAQBvPTCOhtHEkXO92SEKXBmoNEyYTdqzIGFwcpMIEuHoBpO4DcApSukZSmkDwNcA3OM75h4AX7K/fgjA26KYASFkM4A8pfQJSikF8GUA7+3Cua462O5lJJPomdzEgkRpDWVP1Y0mEqp1uQ2lNQBAVRfZJKuN41NFXFysrvZpBILdH71inA17w6bKEhKKLOokQtCNILEVwHnu+0n7scBjKKUGgCUAI/bPdhNCniOE/IAQ8mbu+Mk2rwkAIIQ8QAg5TAg5PDMzs7J30gewC3Ekq6FUN9Bs0jbPsGZiP3bscuzfUWVBYk0xCUtuAuAwCpFyuPr4+Fefw59+5/hqn0YgGIMo9mgzxUu/IgU2HN0IEkGMwL/yhR1zCcAOSunNAH4TwFcIIfmYr2k9SOmDlNJDlNJDY2NjHZz26oBVXI9kNFAKlBrtF/KvP30OH/nSYVRiHAusUblJd+UmFixEG4TVx0KlgcU1WtTJ5NheMWKP3CRSYEPRjSAxCWA79/02AP4ubs4xhBAFwACAeUppnVI6BwCU0mcAnAawzz5+W5vXXJdgnsRINgEgXtV1uWGCUuByIV47CxZMyjGDSjcxW6rj4SMXWh7nPQkmO4lsktVHtWHG3nz0G4xJlG1vottgcpMiEZECG4FuBImnAewlhOwmhGgA7gXwiO+YRwDcb3/9fgCPUUopIWTMNr5BCLkKlkF9hlJ6CUCREPIG27v4MICHu3Cuqw6dk5uAeP2b2I77cqEW63e4clP/6fNDz0zi1792pEUiEHLT2gOlFBXddK6XtQbes+uFdGqYTagyASFEpMBGYMVBwvYYPgbgOwBeAfANSulRQshnCCHvsQ/7AoARQsgpWLISS5O9HcALhJDnYRnav0opnbd/9lEAnwdwChbD+KeVnutagONJZKwgwTOJh49cwP/53RMBz7Fu4rhBoqKvnty0WLHej3/h4Y1rR24SN+WqomE2YTYpyms2SLjXby98Cd1sQpWta1JkN4VD6caLUEofBfCo77FPcV/XAHwg4Hl/B+DvQl7zMIDrunF+awl8dhPgvRH+7tkLOD1dwm/+zD7PcxiTmFrqlEn0P0iwm7nm8xt4T0JjTEJ4EqsKdp30k0nUdBNJVY51LL+B6sW1rJsUimTZn5ZxLa7HIIiK6z7DMa4ducm9ES4sVAIvVPZYx57EKgQJFvT86a283CRLBIpEhNy0ymAJDv3yJKaWarjhP/8LDp+db38wvFJsL2oldLPpbFgSiiyCRAhEkOgzHOM64zWuKaWYXKgGLpwNozNPYjWzm9j7qbUECZdJAJYvIW7K1YUbJPoTrC8sVtAwmzh+uRjreI8n0aMgweSmhCKhITYtgeiK3CQQH8yTGLaZBNt5z5YaqBvNwDxfFjim1oFxXahFBAmVCxKqLJjEKoNdJ3XD8iZkqbedb9j1OFeK1wG4WDMgEaBJe9PnTDcpFNmVm0QxXTAEk+gznEEnqoy0Jjty0+RCBYDFGqwicxfs4l1PTKJFbtJduQmwmYTwJFYVvMzUj+p3xgbmSvFk00JVx0Q+aT23J56En0m03nsCIkj0HbrZhCwRyBJBLqk47cInF9zWCH4Zhi2m04V6rIuYZTdV9d7kl0eBMaMW41rITWsOFS4wVPqwoWCbltlyXCahY/NA0v66R54Ey26SJTQpYPT5flkPEEGiz9BNCtWmuPmk6tDoC1z/HD/tZYtpw2xiPsYNVuV2iP0uqAvyJAyzCaNJfUxCFimHqww+q6kfvgRjA7GZRM3ARD4JWSI9SoH1yk0AenJN1nRzTXU/6BQiSPQZDcOluPmU6uyQmNwEtKaGNoymoxfHyXDib/h+Xpw13XQCGh8kWNDzehKiV85qo9LnIOEwidiehI58UkU2ofTFuAZ6EyQ++c2X8MDfHu766/YLIkj0GTzFzScVh0l45abW9NEtgxbtjuNLVBsmBlIqgP4GCV4S4IMEC3pCblpbqHo8id5fJ6xPWXxPwkA+pSCbUFYsNzWbtOXe8RbT9a7A8+R0CRcX4/mJaxEiSPQZ/IWZS6qOPDO5UAVrnu6/UBtGEzuG0wDiZThVGibGclaKbT8znPgMlCrHhtj78ctNIkisLvrNJBgbWKjoTnO9MOhmE1XdRC6pIpdUUFzhZueL//4q7vhvj3vMel767aXcNFOsr9nWJ3EggkSfoZvUuSDzKWuHZNVIVLB1MAWgVW6qG01sG7SCRDsm0WxSVHUTY3YDwX4yCb5C1sMkbGbUyiTW741zJcArS/ZPbgKA+Uq05MSYQz6pWEFihZ7E/3d4EjW9iSXuGjV8bTkAtwVOt0ApxUyxvmabKMaBCBJ9RsNuKga4xvVcuYGa3sTVY1kAQXJTE9mkgtGs1jZIsFRGl0n0MUjwcpPBB4kQT0KkwK4q+LTXvshNXCCaLbYLEtZibjEJdUXX8SuXCk4BHx8YGyaFInk9CX9W3kpRqBpomM2uv24/IYJEn6EbXrlJNylOTZcAgAsSrXKTpkgYzyXb9m9iNwELEv31JDgm0QjyJITctJbA7277ZVyzBIy5crQv4TCJlLpiT+KbXOt6Xvax2nL45KYuF9RNF2vO67aT2NYqRJDoM/h+MfmUVfD+8sUCAODq8QwAb5BoNikaplVjsGkg2Ta7qdrlIDFTrGOmGN9oBACJeHdkQm5am6hwCQ790MzLDcORVNtVXTPpMmfLTcvNbmo2Kb515CJySeteq/pSs53sJrk3ngR/76zXcb0iSPQZllnGspusG/SVS1aQuGrUYhL8hcp2NpoiYSKfbCs3VWzZgHkSKzWu/9NDz+O3Hno+1rHMuB7NJoLlJi5IiK6bvYVuNvGhLzyJpyOa6VUbptOyvl/G9c4Ry1ubbZPhxKRLK0ioy2YST5+dx8WlGn72Zmv6ccXDJDi5Se1NkJgWQUKgU/CeBNvdvHypgLztOQBeT4LPDJrIJzBXbkTuwNlNMJzRIBGgVF+Z4XdpqdYBk9ChSATDGc2zM3WYhCracvQLs6U6fnhyFj86ORt6TKVhIpdUoMlS34rpNg8kocoEc22KQtmGI29nNzXM5rKY5zePXERKlfGem6wgUfVkN3Fyk9ybFFgPk1inGU4iSPQZfAps3qb6Jy+XsG0oHTj7mZdqNtl9bKIWbXYhpjQZmYSy4qyVYs2IvYAUajryKRVJVUaNu9mC6ySsBn+iV05vwKS/6TbXSkqTkdLkvmTflOsGsgkVI5kEZttsPNzsJtXZTHXKJhpGE4++eAl3XTvhbMAqPk+iJbup20GiJJiEQIfwFtNZQaJhNrFtKBU4+5ldtJoiYWKgfUEduwnSmmxVqnKeRN0w8fr/8j18+4VLsc+3UNVj+xqFqoF8UkFSlbzGdYDclFBEr5xegiURzBQjrhXdQFpTkNHknjOJpj0BL5uQMZrT2jMJ25PIJq1iOqDzduHfPz6NpaqO9960FSmbxYbKTT1KgZ3m7tV+tWTvNkSQ6DN0g/ck3E7tW4dSgbOf+QV2ImcFiaml8F0Y2xGmHSbh3lgzxTouF+o4O1eOda5mk6JYN+IHCZtJpFTZ50kEyE0BAVGge2ByTRSTqHBMotdSCGsmmEkoGMkk2lZdF2sGsgnFboSpOo91giPnF6FIBG/aO4qUZl17/uwm1Zfd1G0JdKZUd4pkayJICMSBdWF65SYAXrkpgEmw7CYgmkm4cpOCjI9JsJ45cRdmtnOr6CaaMXb8harVayepyr5iumC5CbBaiAt0H2xBnY7Ihqs2TKRVGWlN6bncxK6lbFLBSFZr27+pUNOdTRRjEsUO/bVizUAuqUCVJaQ16zX8cpPWUkzXZeO6UMdmWybe0EyCEHI3IeQ4IeQUIeQTAT9PEEK+bv/8SULILvvxnyGEPEMIedH+/63cc75vv+YR+994N851tcEb1wlFci7SbUOpwN0Mb1wPpVVoshRPblJlZBOyhwUwHViPeSOw3Sil3uK48OOtXjtJVfbor2G9m4Du35QCFphcM1uqhwb4SsNE2vEkeruAsc1KNqFgNJvAXDm67X2xpjsMYrmeRKluIGs/V5YINEVysv/MJkWTwpGbeulJ7LAzujasJ0EIkQF8DsA7ABwEcB8h5KDvsI8AWKCU7gHwZwD+xH58FsC7KaXXA7gfwN/6nvdBSulN9r/plZ7rWgC/eyGEODfAtqEUZIlAlb2zn9lOO6FIIIRgPJ+IZhI6Z1xrXuOapR3qMW+EpQ4H0XuZRJD5HiA3iQynnoClkBpNGtoCwzKuFaQ1uecLGNusZDQFIxkNNb2JckRgYs39ADdIdOpJFGs6cgmXrac5WY1tlJjc5Eq93bse64aJxYqOncNW/dNGzm66DcApSukZSmkDwNcA3OM75h4AX7K/fgjA2wghhFL6HKX0ov34UQBJQkiiC+e0ZsHXSQCu5LRtyNpt+CuR+ToJANiUT0Y2+as0rJGPCUVqMa5ZkIi7e+cb9lViZEm52U1Si9xECBwGBSBQWhPoHvjPLkhyMswmGmYTac2akNjrynwnSCQUjNg1PFG+RLHOMwnmSXQuN2U53y+ttgYJfugQ0N3rkUlqG55JANgK4Dz3/aT9WOAxlFIDwBKAEd8x7wPwHKWUv3L+2paaPkkICRzASwh5gBBymBByeGZmZiXvoy/g23IAdgOzhOJUvmr2GEUGf0sLq6Au2ozMaAoIIZZx3Wj1JGLLTdX4w4vqhoma3rSzm1o9CcaEGIJM+uXg+fOLOD5VXNFrrFXMlxt4mGsp0Ql4aWY6IMOJGclWkFB6vstlXVxzXD1QlC/BMuUAzpNYhtyUS7hBIqXJzvvWTUvqUuw2IYQQaLLUVbmJparvFEECQYu3X2yMPIYQci0sCepXuJ9/0Jah3mz/+1DQL6eUPkgpPUQpPTQ2NtbRia8GGlxGBQAMZTRnpwG0tqvwMwlWdR2m57LcdwCt2U32zi3ubsnDJNosInyvnZQqQzep06vGP98a6A6TuLRUxQc//yQ+849Hl/0aaxl/8+Oz+PWvHQmcRtiuvqRQ1Z3dcVCGE19Pk+YWz16BZxKjcZgE50loioSEInXc5I8Z1wx8MDQcucnrk/UiSGy3VYKNbFxPAtjOfb8NwMWwYwghCoABAPP299sA/AOAD1NKT7MnUEov2P8XAXwFlqy17sF7EgDwe++8Bv/t/Tc63/uH8fj7Hm0aSKDSMEP76zMzErB2bbpJnddwjet4tQmFDjyJItdGIWn7Daygzj/fGli5J0Epxe/+w0so1Y11PdAlCqynlz9IHL24hP2f/Gecn68EPQ2A9XnsGrUWp6DiS76eph/GtRskZIzYTCKsVoJSioJvgc8lVU+X4TjgjWsAnqJBtvlSJX+rmO79HRiDm8gnWyTY9YRuBImnAewlhOwmhGgA7gXwiO+YR2AZ0wDwfgCPUUopIWQQwLcB/A6l9N/ZwYQQhRAyan+tAngXgJe6cK6rCpZRwctNe8ZzOLgl73yfUGRvdhOTm+xFdSRj7cLmQ6h6xTYjASBjBwtmXjueRMwbgb8p23kSLKAw4xpwZ0rUjaanTTiwcrnpkecv4rFj0xjLJTC1FM6s1jNYT6+lqvezPnG5iIbRxFE7iAShUNMxnksil1RCgoT12aZUBWlVQcNowuxhYSPrIZZNKBi2+0WFVV1XdRNmk3pSxHNJpWMmUapZFd4MKY8nYb1XntX7pd6VYqZo1UiMZLW+pBn3CisOErbH8DEA3wHwCoBvUEqPEkI+Qwh5j33YFwCMEEJOAfhNACxN9mMA9gD4pC/VNQHgO4SQFwAcAXABwF+t9FxXG05GhRz+Z/fPfm74DLahjHXRL4RlrOiGwyQyth7rny28HCbh9yR+dHIWf/79U+6xrNdOKihItMpN2gqySeZKdfz+t17GjdsH8Utv2o2qbnr8kysBS1UdFxatkbYLZa9hO29/z34eBCa1jOcSgZ5EtcF7EqwauXd/w3LdSqhIqTISioxcUgllEuyz9DKJzgYP1XQTDbPpk5tcxmQE3IuaInU1JXu6WMdwWoMqS3aAWp9JGkr7Q9qDUvoogEd9j32K+7oG4AMBz/tDAH8Y8rK3dOPc1hIciisHevAAAuQmh0lYN/JQ2tqFhQWJSsN0jD6nnUHdQMNwp3LF3S0Vajoymoxyw2zJfvn7Zyfxjy9cwgNvvgqKLDk3diCT0APkJseT6JxJ/NGjr6BY0/Ff33cDTk5bpvVUoYaBtNrmmV48fnwan3r4JfzLb9zh+DhrBYxFAK2f9aL9/eRClNxkpSOP55KB2U1+uYk9xnyAbqNUN5BJKE7ywlg2EdoJtsg192PIJjprF16qtwYaXlZj96IieT2JbqZkzxTrTsv+pCr1ZbBTLyAqrnuIz3zrZTx+3C3vYPUJmhLBJMJSYBmTYEGiHLyrqjZMp08NzyR4XTt2CmzVcKq8/Zp1oWZN3JpcqNrfMyahIOmb8hXoSSyzDQKlFN89ehk/e/M27N+Uc5oeXloK31WH4VvPX8T5+WqsueH9Bh8k+HoVwA0aFxbC33OhajOJfCLQuK5wxnUm0drXqNso1Q1n0wJYEkzYTImCM5XOzyQ6CBKswjvhZRJVX3aT5pebusgk+CDRjwyyXkEEiR6hppv44r+/in995bLzmKODRslNPvOsrpueGoM4TMIvN5XqhmfXFj8FVsdoNgFVJi1Mgt3Ip2dKzrGAtftjO9MouWm5xUuXlmoo1g1ct9XyceK0KgkCpRQ/OT0HoHUR7geOnF/E48fC60NfuVTAcEaDLJGWz3qhjdykm01UdRP5lOrITX7Phu1q05qClMpaVvRWbsrwQSITziQKXKYcQzbR2QhTN5GCL6YLyG7i5aYepMCyIJFSe1+w2CuIINEjnLMzT3gdMo4n0VIn4asxyCWVwIWDgTeusw6TMJ3016G02pHclE+ptunmYxJVX5Co6ZAlgrQmO3JTNcq4VpcnN52wZxXvncgBAMZzjEl0FiTOzVdw0X7OagSJ//LoK/i9b4bnYrx8qYCDm/MYTKlYqIQwiZAgwWeajeeSqOnNlmy4SoAn0cudbskfJLLhnWDdDYeXSRQ68CRYnyeeSbCFmk17BPxyk9y17CZKqTdI9KGJYq8ggkSP8Oqs1WmV1yGX5UkY3pRZSSKBCwdDtcEb1yy7yXAySbYMpjpiEvmkiowmt+zi2EJ0erpsH2sVPxFCHLnLkZsCPYnlyU1sHvg+O0hoioTRbKLt7G8/nrBZhHXu/Q0ShtnEi5NLmCrUAjOKDLOJE5dLuGZzDoNp1fEgGNhnv1jRA3fXvKY/nrcWKb8v4a+TAFYmNz3z2gLu+dy/h1ZuW7MkXDY5kk1godIInPvMz5JgYNlNcbPYSrVWT4K9z5phwgiTm7rEJJaqOhpm05kQKZiEQAtes9tx+1sTA/As+n60pMAaTU+LbQAYTKtYCCmwquiu3JT1yE3W8ZsHUh0wCat/TjrRmr7XIjfZrAOAWycRITcpEoFEOpebTlwuYjSrOWmU1nuKblUShCfOzDnn2W8mceJyyUnzDEpPPTNbRsNo4prNeQylNSz6mUS54aQ3B/kSfHYQW6T8GU58I8hUF4LEC5OLeP78Io6cXwz8ud+TGMtqoBSBm50gqSiXVEApIvs9Bb9Ga5CoNMxAVs+P1G02Kf7631/FUshmrB3Y5zpue2bpPtSi9AoiSPQIr85acpOnNbERw5Pwp8D6mARgjSYNkpvqRhOUwlNxDdhMolRHWpMxkFJjpcAaZhOlumExCd+Eu2aTOjtY3pNgOz8WEDxyk49JEEKWRe9PXC5h73jO89hEPtkRk6CU4onTc/ArMIcAACAASURBVLhjn1Wh3+8gwS+kFwMMd2ZaH9yStzYE3EJFKcVCpYFrtwwAAC4stmY4Fbl0ZMYk/MGo0jChyRIUTxvt5XsS7LN+fjI4SJTrpk9usquuy61BslDTocrECeKAGzDiZjjxXWcZmAxb5YKEP7uJsf0XLizh97/1Mr79YvwBXTzY35sF6aQmb+hiOoEAMCbBXxiNgFYAfgRVXPv1/MGA3SXg3R0CVjDS7HYGs6U6RrMJz24pCuwmy6csuYmXEUoNA5RazQYXKjrmSnWHdQBukKpHeBKAFRA7ofeUUpyaLmHfRNbzeKdM4sxsGdPFOu7YN46EIvVdbnru3ILz9cUAX+HliwVosoSrx7L2Z+1uCKq6ibrRxHVbrSAxGcQkuOygMduzaZWbDOdz6obcxBjz8zGZxIhTUNfAYqWBJ8/MOa/BWnLwvb7c/k3xPit2nKfiWuWZRLDcxFj80YtLAICpZWTNAW4rFN64jvr7lupGKAtbbYggAeAnZ+bwx/90DA8fuYATl4uxNfsovDYXwCRieRJyiyfhl2qG01pgPx93Kp17Y7BOsLOlOkayGjSZxHp/bt2DgrSmeGg+o/I37xgEAJyeKXuYRItxHdC7yXqv8QIWw8WlGkp1wzGtGTYNJLFY0WPv1Jgf8carRzCQUleFSdyycwgAcCmgpcjLlwrYM56FKksYSqse1shYxd6JLDRZCpabOE0/n1SQUKRAuSntCxIrMVbZc1+YXGr5GaW0NbvJ3mH/2leexU2f+S5+7sGf4JMPW0Y+S9/l4cyUsDcrNd2M/LyLdcPu+eRed3zRYJDcxDMJVs2+3PRoV25iKbCWJxHmqXztqXN4/1/8eE1WZXelmG694+jFAr7wozPO7mI4o+H7v3WnxzjrBDXddGQE3qyK50lIMJtWczzFTsnz11UMZlQsVnRQSj27LT73nSFjDx6aLTawYyQNVZbiBQlOssgkZM/Fy3beN+8YxD+9NIXTMyUUa4YbJGLUSVjvVe4oSLDMpn3+IJFnY11r2DWaafs6T5yZw6Z8ErtG0n0PEsWajlMzJfzG2/bh2KVCiNxUxJ37LSlsMG3NXqjpJpKq7HhRwxkNWwaTmAxgInw6MptB4q+VqOgmxyRap7Z1CtYg8NJSDdOFmqPFA9bnbzSph0lsH07hLfvHkNJkXLd1AGdmynjomUncd9sOpxCQh3/w0C9/+TAUieCv/6fglm6lmrcDLADPCFNHbgpJgT16Ycl5P8vBTKmOhCI555BUZVBq/S2SauuGabbUgNGkmCs1kB5eW8vy2jqbVcJH3rQbH3rDTpyZLeGbz13EX/7gNM7NVRxK3ynOz1dAKZBLKJ7dTty2HIB1MSmyZJu+3uOH0hoaZtNqC87dCHxaI0NGU1Cqm5gr1/G6nUOxMzj4hcbvSbAb9cCmPBKKhNPTJRRqunMjK7IEVSao6SYM01ogwplE/IXppBMkvHITq5W4FCNIUErx5Jk5vHnvGAghyPc5SLwwuQRKgZt2DGLzYKqFSUwXa5gt1XHNZqsOZNCuIl+s6Ng0IDusYiitYdtQOpBJsM+HSS1BVddVjkn4p7YBVsBtGE1Ph+Io8PObn59cws8cdINEkD+QUGTPAl+uG/jRyVl86uGXoMpSAJNwPYlXLhXww5OzGEyrLRsl/m+Q9b2G17hm/iDXvl61PDLDbOKY3X6+0/obhulCDWO5hHNuPFsLChJMzp0vN7B9ON7fvF8QcpMNTZFwYFMeb792AkBwD/64YOmv12zOe3ZnjRjGNWMZbCFvBOzCh+2COr/k5DRt44JENmHll8+XGxjLWn1kjCZtO7Oar6D2exIsgAymVewezeD45SIqDdNT/JRULHrN6HuYJ9FJCuyJyyWM5RIYTGuexzspqDs5XcJsqYE3XmWNM+k3k2C6803bBrF5INnCJF65ZC1O12y22JK/eJLJTcMZFVsHU6GeRDZh1dMACOzfVGkYSKu+amTuWv3kwy/h17/+XOz3VdVNbB+2piv6fQm+TXgYMgkFv/vOa3D0YgFHzi+2MAnek/jyE68BsAJnWK1Fqd4qWTkLtc5lN0mtTOL0TBl1o4nhjLZsJjHN1UgArh8SlgbLAmnYFMHVhAgSPjCaHDXYpx2YH3HN5pxHh3TkJiXCk3CKzMKlGn53ycNt2ubeHJmEggsLVTQpMJpLxB74zvdiSmuKk7IJ8MaoiqvHs87Cxxc/JewRpkHzrZ1jOpSbTl4utrAIAFxrjvY39PftNilvvHp1gsRz5xZw1VgGA2lrkfe3OXcym3xMwgkS9qI4mNawdSiF2VK9RZv3z1GwgkQrk+A3E2lV9rDFs7PlwJ5PYag0TAynNeybyLVkOLlMIro/1rtu2OwE7zBP4sJiFd987oIzyIfVzfhhdYD1y02t2U18EommSGhSN0Przv1jKNaMZU3tu7RUw5bBFPe7o5MDGPsLSm1fbYgg4YOTV76CIPHqXBmDaRUTA0lHhwRiyk2+FtpBxvVQJrg1R5DclE0ozm51NJtwmEo7X4L3JNjNxnZBbrGTgqvHsp6BQwwpTULdzsSx3lfrAqHJ8eWmZpPi5HRr+itgBcJcUmnLJJaqOv7yB2fw+t3DDqXvZ5CglOLI+UXctN0y/DcPWIs8/zc4dqmATfmkw5YGU9b/LF+ffeaDKRXbhqxFyJ8h5df0x/NJFGuGJ5jwxjVgVwTbchOlFBcXqx39Xaq2Z3LT9gE8f37RY9Cy4BPFJAArLfr377kWikQcY5vBmrYIfPWpc6jqJn7vnQcBWMwwCAVuaBFDWuWN6wC5yb73nju3iKQq4aeuHgXQuXndbFJcWKg6nw8ArsA0jElYf+ughJTVhggSPmiKhOGMhssrkJtemytj10jGpZj+4ettiukAN7AEGddh/ZucKlrVGyTY/cpSYNnrRqFQ1UGI5aukWQM4e0fF5KZcUsXVY64HwC9MTG7yD03yvFc1fnbThcUqKg2zxbRm2DyQbNvk77P/ehILlQY+9e6D7jmnVBRrRk9nKTBMLlQxW2rgZhYkBl3DneHYVNGRmgC+NbxbZZ1PKlBkCVvtnapfcvJnBzHZg9/4VHxMIpNwW68UagbKDROluhE7069mF3HeuG0QhZqBs3Nu/UYpoEVGGPZN5PDNX/tp/MrtV3kelySCrKZgttTAoZ1D+A/XjCOjyTgdxiTq4cZ1JUxucoLEAg5symOL/flc7lByminV0TCbztx6oH1yAAukIkisE4znEitiEmdnK9g1knaNMp21J47X4A9w21UEG9f2whHiSaR9Nz/DqO1JAO1nShRsui5JBBnNrdwGrPTCpGrVYFw95so/Hk/CnnPtMIkgT6KD1sysJXiQ3ATYBXURn9mp6RK+9OOzuPfW7U4hGgBntngnswqWC8eP2G6lv24ZYEzAWoQsPbyEA5vdIVT+DcF8ueEwya32TtXfw6lY1z2fxXiuteq6qvuYBJfHzzOTuDUkLOjcsM0KgC9wkhM/cCgOrts64LxHHsyI/vBP7QIhBFePZ8PlpnqrcZ1QJEjElZtkiUCSvHUSAHD8chHXbsm7WXMdMgnWwt3DJDTrtdt6EiJIrA+M55PLNq5Z+uvOkYxbL8CYhBEjBVb1yk1BTGIgpYIQYN7nSbjD7fk6CW+/HPZabeUmru7BX2zF/+wqnkmkvIVLXk8iKLspfsX1icvWYuCvkWDYPJCMLHz6o2+/jJQq43+7a7/ncRYk+jG06Mj5RSQUCQdspsB2qowBnZktQTcpDmxy32NSlZFUJaegbqHScALHpnwSskRaMpz8TII1QeR9CSu7Kdi45oNEXMmJZe3sm8giqUqewrA4xnUcWMWBCdx97SYAwJ6x4CBBKW3xZQA4fcUqDat3k79eiV2jlFqBis+a6wSM2W3ngoS7FgRfZ0yyFUFinWBiBUxicsFKf909mnFuwha5Kcq4Vtob14osIZ9sbfxWbVhtxfl2BuzG1GQJ+aTi3BjtZB6+F5PbTdZmEtwNmNYUR/bg5aaEKrWXmzoopjtxuYiJfMJZ1P3YlE9iplgPbBj3gxMzePz4DD7+tr0Y9Wnd7PX64UscnyriwOa8w+Y220yCLULHnMymvOd5gym3wn6xojtMUpElbMonW5lEiyfB5Cbr9zSbFFXd9MiS/HjNi9yiuBjz78LkJkWWcN2WAU+GU7eCxMffthf/9X03OBudq8ezmCrUWlhgTbdGsfKjSxlSdkfjhtn0SE2Ad87LtVvySGsK8jG8Lj9YkNg62Co3hTMJ4UmsK4znE5gp1ZelU7OeTTtH0i1pb8s1roOGFA2lWzvBVuyBQ3zeOLsxR7MaCCHxjWu7qysApFmQsBcRPoAALpuIlJtW6EmcvFwK9SMAYNNACk0KpyU6j2+/cBFDaRX3/9Sulp/1Ikj84wsX8dKF1srjcsPwZIClNBlDadVZ5F+ZKkCVCXb7aj34/k3zZZdJAJbkxE+oo5Si4NtFD6c1aLKES/ZiVzNaExz4qW3LYRLs2gOAG7cP4ujFgnONMSklo0VnN7XDu27YgrccGHe+3ztuSY+nZ8qe45w24cnWoGQxJstr8bfHYfeGLBHnWts8kFoGk6hgNKt5PB/Xn2y93g2z6RSeihTYdYKJfBJmky4rqrOeTbtHM1zam3WTME9CkcKZBG8sG6a1IwqSaoYyWoAn4dWZAZcFjNq6dGzjuqY7Cyi7uZm5Zi1CbkDYN5GDpkieRSDlBAmbSQQUEFkdb+PJTWdny7gqolBu04D1/oJu6DMzZecc/ehFkPidv38Rf/Pjsy2P1/TWTLUtgylcshflY5eK2DOea9lEDHH9mxYrDY9ev20o5ZGbWKoyH7AliWDHSBpn7fqdoCy4tC9IsEs0jidBqc1M7N3yjdsHUTeaOG4XpJVqBlKq7Klu7gb22EGCFVky8Nl3frD3GSw3Wee3dzzryEMTA8llMYmtQ96COP9awIPdV4RcwUyCEHI3IeQ4IeQUIeQTAT9PEEK+bv/8SULILu5nv2M/fpwQ8va4r9lLMKNvOdWWr86WMZBSMZjWWtLedNPq6BpUIcrAT2xzCtECmURrJ1i+aRuDyySs96TGZhIuW2BMgl3gxaruuQE/eufV+PIv3uZ5X0lVilEnEY9J1A0TxbrhKU7yY1Pekm6CMlHOzJY93gkP5qN0K0gUajqKNSNQVqgbpkcKBLw71eNTRVyzqZUtDWVULFZ11A0T5YbpyE0AsG0whalCzfk8+TbhPHaPZnBmxtu+PuXxJNypbZcWa7jKTkgIaiTZ+r6s382ud5a99azdyLDcMFYsNQVhx3Aamizh1IzXlwgaXcqQ0twiT8UnN7Fr9OAWV+7blO98VsmkL/0ViE6BZcxncz6JpaoeKJlGwWxSHJsq9Cz5YsVBghAiA/gcgHcAOAjgPkLIQd9hHwGwQCndA+DPAPyJ/dyDAO4FcC2AuwH8OSFEjvmaPQMrqAvq9d8Or81VnNYQfIUnYBnXUc39AK6YTm86u/2gHbA1jKZVbuKraAF38NBoVvO8Vnsm4fZiymruhDvrZ94c9NFsAm+wi6AYUqqMmtFGblJkGE3aVtZj4zqDMl4YNoeYjIuVBubLDVw1GpwV1SmTOD9fiaxWZ7v6IIZUD2QSlqewUG5gqlBzTG3vOVpMgn3e/N9h65Als7GFjB84xOOqsQxem6vAbNJQJtEwLfZ6YbHqmOdx/i5OzzA7AG4bSmE8l8Czr1lBolQ32xbSLQeKLGHXaLolDTaoDQgDM+gNk7bcV+x7Pvtt00AKM6V67FRgp0Zi0BskNEWCIpHAFFh2vtuH06A0vg/EMF9u4O7//kP8w3MXOnpeXHSDSdwG4BSl9AyltAHgawDu8R1zD4Av2V8/BOBtxNp23gPga5TSOqX0VQCn7NeL85o9w0qZxC67GtRfZRmkg/rBexJRhWhBnWCrutnCJLIhTCKq4tqZJeFr/c0MSL4teBiSqnUzRspNaryAxWYOjEQEicG0Ck2RWj4zpleHMYmUKkOVSazF8HKhhrf86ffxjcPnQ49xgkTAe7Ka9LUyiWLNwDP2gnpgU77leUP2hoB93rwnwXLxz9u+BN8mnMdVoxk0zCYuLlYD27ekOUlxqlDDzpE0Mpoc6+9S9WXVEULwuh1DeIYxiXpvmARgSU7+DKci1xHAj5SqOEOH/Bu23aMZHNiUw1vs5oqAlRBBafwNo1sjkWr5Wdh0OnZf7bALPDutumaKwnDE/bESdCNIbAXA3zWT9mOBx1BKDQBLAEYinhvnNQEAhJAHCCGHCSGHZ2ZmVvA2XDjFRx0yCd1sWumvw94gwWh8w6SRpjXglZuYVBNoXGc0VH3tkoM8Cb/clIjBJJxZEvZNpikSNFlCuWH9vobRbNsh12qW5hpyYXKT9V6jfQnGJIYz4XITIcQuqPMGiTO2FHHVWDCTIITErro+enEJRpPie69Mhx7DTOggWSGoAyhLg33MbhcSxCSG0hqMJsV5e276ICc3MfOW9XwqBFS/A8DuUWbyltz2Ldy5sGv17FwZZpNi80AqdG6JHyytM8lde7fsHML5+Sqmi7XAFhndwp7xHM7NVzx/76CpdAysZbceIDeNZBP459+43XOthDHUMLg1Eq1N+sLmXLPzZa1G+H5Ur86W8fOff7LFd+HBNg/D6bUbJIL0Ez8fDzum08dbH6T0QUrpIUrpobGxsaBDOkZCsbJOOmUS8+UGKHXlqqCK66gaCfa7AeZJhKePst0kfxMHBYltQym847pNuH2f1WIgTjGd07eJW2jSdrvwKFOQB3vvzPgMk5uA9um4jEkMZ6ID00S+ter6zGwZqkw8Oet+5FNqLIP2+JQVcJ44PRsaZN0gEcwk/H8H1t/n8WPTGMloTlsYHgN2UGCNI/kd43g+iU35pFO85nbvbfUk2Gu4clPraE/W5mLrYCp2h1yWscOn1L5up+1LvLbYMnCom9gznkWTWsGNoZ3cxNpytGP1gHVNAfFVBZb+GsgktGAmwctNgJdJ/Pj0LH50ahb3/dWTODUdHCgchrmGmcQkgO3c99sAXAw7hhCiABgAMB/x3Div2VNM5JMdM4lZO/2S6f+qbOmQfApsECvgocoEhFiadtQunJmXvORkGdf+KlMZf/Hzt2CP3fMoTjGd07eJW2isluOGp6dTFJissugEieBW4QDaVl07O6UIJgFYu+pjl4oej+PMTAk7htORmTVxmQSbZ1FumJ7pcjyY3ORnEqxlup9J8DvV/ZtygUkNbEPAgsSQb8d4w7YBvGgP+wmaDw1Y12QuqVhBQm+dO5KyvSw2jnbLYAoDKaUleDabtOXaceUm9/Wu3TIATZbw3LmFnhnXgFVQB1gp0gz+Vuk8knYxnW42oUZkGTIEMYkXJhdDZ187NRJhclOQJ1Hzyk08kzg3X3HWhHsfDA4U8+W1Lzc9DWAvIWQ3IUSDZUQ/4jvmEQD321+/H8Bj1OoA9giAe+3sp90A9gJ4KuZr9hRjuYRTfBQXsyXrw+ILtvj88yAd1A9r9rPkyW4KNq4Zk3AvKMu4jjYI2e+Pkpuc3SgXCDIJGZW6GUnlebDFcMnuARX0vrXYclMDhCC0kI7hlp1DKNYNp4UHYC2su0NMa4a4QeLYVBG37ByCLBH828lgaZMNAar53lOYgT+RT4LFhSA/AnA3BGfsIMHLTYAVJM7MlrFU1Tmm5z2GEIKrRjN4dbbsyEN+4xpwu6puGUxaRXxVrz7+hR+9irf/2b95HmMeBx8Ak6qM67bm8ey5hZ56EleNZUCItxtsyW4bEyTtMuPauhfbL39+r2u6WMN//PMfB6Y4A5bcNJLRPCyNYTlM4txcBduH0/jqL78BgBUoZn21QG5n4OUNSWuHFQcJ22P4GIDvAHgFwDcopUcJIZ8hhLzHPuwLAEYIIacA/CaAT9jPPQrgGwBeBvDPAH6NUmqGveZKz7UT+JnEdKGG+x78SWQTuTn7w+M7WLJ6AcCaJxHnwmQttKNaWrBdA19842//HAQthnFdCMiQsUaYGp5hRFFIcUEioQSn/fL+SxTm7AIyuc3Oj40EZSaw2aQ4O1fxNCEMwkBKdd5zGAyzidPTJRzaOYTX7RjEv52YDTzOzW7yvid2DfiZhCpLmLDbZgT5EYC7ITgzU0ZGk1uuB9Yv6eiFJRRqOhSJtBjkgJsGGzicys4+Oj1dQi6hIJdUA4Pn0YtLeHWu7MnwYu8t5Xtvr9sxhOcnl+w+YN3PbgKsv+f2obQnDbZY0wOrrQHrPRt2hlccuYkQgk35pJM99q3nL8Fs0pb0c4ag9FeGUCZhB4mhtMX2/Exix3Aae8az+Oy9N2G2VMdz57yt2OcrDWQTSuA60Q10pU6CUvoopXQfpfRqSukf2Y99ilL6iP11jVL6AUrpHkrpbZTSM9xz/8h+3n5K6T9FvWY/MZ5LYKZYd26G775yGU+cmcMPTwYvDoArN41kXdqXbmES7f/kms0k2A47rOIacLuDUkpR0U3nZo96baAdk2CehHdWdqVhemZJRIEtUksVPfTidWdnRDOJ+XIjFpXeMZzGaFZzgsSFhSoaRjM0s4khDpM4O1dBw2xi/6Yc3rx3DC9dXHI2BQw13XSuAb/cxAJh0OLNusFeE8Ik2A5xtlRvGbgEANfbExSfn1yyWnKk1MCgvHs066TbAsFy02vzFccnGUi3/l0uF+qgFChxRWFBchNgBe2GYaVyhy3a3cDe8SxO+eSmMM+MybGFmh5LbgKsoVYsSHzTTjNljM0PK0gET5bj1wIerNhQlgiGM279E6UU5+YqTiIMS60PYhK9kpoAUXEdiol8EkaTOjv1p16dBwCnijQIc6UGNG6uLWCngjpMor1xDbhjPRsRNQaO3GTf8A27OjuI5vKIU0wX5Duk7el07uyI6N+T8DGJwGM68CTiZG6w1EuWn396NjqziWHANq6j6h/4+dq37xsDpcCPTnk3DEy3Hs8lUDOCmURQwNwykIJEgL0hHW4Huc8haDEYymjYMZzGixcWW5r78WDB8uVLRWtkKXctsgXebFIn42ogpTrztRlYC/0St0gGzVYHgNfZzA5A283LSnBgcw6nZ0rOZiOoAywDe5+FqhFrwwZYabBThRpOTRfxot1uJWgQUdAcCR5JTlXgUW645zuccVPbFys6inXDkaHY5tOfjjtf0XtmWgMiSISCr5Ww5iJbQeJERCraTKmO0Yzm2cXx3TWtOon2u5eEwyTCd5+aIiGbUJwgFjRLIgjxmITlI2R9E+54uakdk2DnsVhtBLYJB+LLTXGZBGDtXs/OVTBbqjsVxlHtPABrMWz6dsd+HJ8qQiJWNs31WwcwmFZbWCWTmq4ay6BhNH2STPhn+Z6btuAjb9odOPsYsIrG2MIfpjtfv20Az59famnux4NlOL18cQlpX48vngU4TMLpkOuyictO0R7HJBrBUtpEPuk0f+xVdhMAHNw8AKNJHfO6GJFy6wSJmh5LbgLsDsOFGv7huQuQiMVYSwFBIqpGgv3uICZRrLmzL/j6p9fmWR8463NLKDIGUmoLk5gv1zHcIz8CEEEiFCyNdbpYx+RCFVOFGhKK5AxID8JcqeH0SGLgzaq4cpPV04iruJaDFw++6jpIZw4C6xsVzSSsi5bvtZ/WXONaIu2btSU9TCJEblI6kJuy8YMEADz72gLOzJQwkFLbBhi2qIZlrADW5mCX3f5dlgh+es8ofnhyxjOB7cKidVOzGRt88HO74bb+Ld5+7Sb87jujGwqwjCZ/ZhPDjdsGcGGxitfmKqFMggWJi0u1ll1/KiJIMMmpVLeGEVlfu3+rWojcBLhsolfGNeC20XjZHv1aCmgTzsA2L5Qittw0kU+iYTTx1afO4017x7BrNBMYJKJqJNjvDjOug5jEOTtIsKwnwMpS8zOJhbJgEqsCZ1BLoYYnbanp3TduwUyxHtqEa65cb6kKTqlKR8V0AOuOGj2wB4BHvwyj/H4QQqApktNsMAh83yaGbMJNgQ3TvHmwHbPV1C6ESajtmUTTNgnjFgpdt9VKvXzm3ALOzJTt7Jfoc83HaM1xfKro6UJ7x94xXC7UcZxjlhcWrMZ4u+ydHx/8nHTmkM+yHRiDCAt412+1zOszs+VQJpFJKJiw24b7F3RepmRykzNL3f678LUCBZ/cpEgk8Nq+ZYd1Xr1kEjuHrQFfL1+0g0TdiDCu3fOILTfZabDz5Qb+481bkE3IIUEivEYCsIoNw1Jg2WAvFiQsP8JiwtuH3dcbyyUCmET8+2M5EEEiBG4P/jqeenUOg2kV77phM4BwX2K22GiZV+BnEvE9Cc64DnnOYNrtBFsNKJAKgyZL0XJTgGSR1hTUjSYWKnrb9FfAK3utxJNYqupo0vg54E7q5WsLdvprtNQEBMsqPGq6ibNzZezjmu+92S5O/LcTbirs5GIVE/mks2vmC+rYZxkmKbWDM/M6RFa4bmveSaWN+nxYDyt/PY0sEefzYBPzHCZRaQ0SvCcR1A6G4W3XTGDfRBb7AxoXdguSRHDN5rwTJKzeYmHGtXuecaRfwA0SKVXGXQc3IZtQAj2JqBoJAEiritMfi4efSdSNJioNE+fmKxjLJTz39Gg24aTaA9Z9X9VNwSRWAwlFxmBaxeViDU+9Oo9bdw07w2CCfAlKqcUkfEEirfo8iTZ1Eux31w1XbgrbffIzJVi/mnZyE2DVLETKTdXW3kzMeLy8VGub/gp4F8OVyE0sHXAkptwEWJLT8+eXMFWoecarhqFdk79T0yU0KbCfYxKbB1LYN5H1pMJeWKhi62CKY1EBTCKmDu4Hy2YLk5tySdXxXqIKHXePeZtP8mCPhclN/CAu3pOo+QYY8dg+nMa//K93OK/ZKxzcnMfLlwpoNqk137qN3ASgpS1HGNgY07uv24RMQkEmoXiCJENUjQTgjjD1JzXw87jZYj9fbjjprzzG7KxLBqYkRPU1WylEkIjARC6JFy8UcHaukGVhdwAAIABJREFUgtfvHsZ4LoHBtBroSxSqVqn/qG8xS9ltAADWBTYek2hwxnUYkxiymcRCuYH//K2jSGtyrEVRU6KZxFJVbzGm2YU/VajFYhKeIBFmXMeQm9hNELY4BuGWnUNOHUg70xpw216EBQm2KfDvhu/YN4anXp13Pt8Li1VsHUo5wY8vqFspk3A8iYjF4Ea7XiKaSUQFCQWEuDvnwZSdQRcgN/FtqSsx6nN6jYNb8ijVDRy/XASl4fIW/77bdT9g2DyQxEfvvBq/9pY9AGzptWF4/CggukYCcNmbf6ZEiSs2ZIv9QqXhSX9lGM0mUKobzsaz1y05ABEkIjGeTzhjGG/dNQxCrIlVQUxitsxacniZRNKe9QzYnkSMC1PjUmBliYS2lBhKayjWDfzCXz+Fs3MVfP7Dh5wbPAqqLLVNgfVXNzMmMRWbSbjn3FZuiggSc6XOWw68boebetku/RVozySOXy5CkyWnuy/DHfvG0TCb+MmZOZhNiqmlmodJ8DJa3cluWt5iys5xKCKL5YZtVr1E1OfD5LegnX9KkzGeSzgbmVzSChrs7zJVqCGjyZAIPJp8tRHOJPqFgzbLf/LMHIDw7Ds+SMRh9YDl4/323QecIUfZhAJK0ZKpdHGxGsmYnJkS3HQ6SqnVADHpZRKXlmq4VKg56a8MrPko8yV63ZIDEEEiEmyAfFqTca2dQbF/IocTU8WWXcRssbWQjj2X6ZDxPQkru6lutDaE4zFkN7x76WIB/9d9N+On9ozGel+aIqEeKTe1ehLMWGuYzbbprwDs4UrW12FyE/tbRAWJ+WXITeP5JLYPp0CI21kzChnNylgKZRJTRVw9nm0J1od2DSGlyvjB8RlcLtRgNCm2DqWcQOCRmyJmfcdBO7kJAK63mURU+xIWNIOYREaTnbnbgKX15xJu/6bpQh0TA5bnUozpSfQL+zflIBHgqbNWkklYnQR/nnHlJj/Yrt/vSyxVjcBiRwb2N6/o7vPqhtXTK+tjEi9OLoHS1uuXNYBkY3qXw7Q7hQgSEWDm9S07h5wFYv+mHIp1wzMsHnC18xbjmptzHduTUN06iShKvHPE6lvzpx+4AXdduyn2+9JkCXrIwmyYTZQbZstCk+aKodoV0gHW7ou997CF0e1TFe5JLPcmeNOeMeyfyMXaubdrF358qoj9AYVuSVXGG68ewQ9OzDjdXz2eBPc3DmvLERc3bB/EzpE0dkQEvZu3D+IP33sd7rp2IvSYbUMpKBJpMa4B4KN37sHH37bH89ggNzr1cqGGiVwS+aTaUiex2kwiqVpSKyt6zYXKTe7jceUmP5icV+SCBKXUzgoMvzf8XaEBl5HlfEziiK1g+D0Jtr4wX6IfTKJ3eWlXACZsavf63cPOY0yXPjFVdAqFgOCWHIB3pkT8OglXboraed6xbwxHPnVX28Z3fmhKuNxUCKmoznA3Vxy5CXA7bkalfWqKFJndNFdqIKPJHS+un373wVijURms/k2tZmShpuPiUs2T2cTjjn1jeOzYNH58ypI5tg2lHHmRZxJRY1zj4HU7hvCD33pL5DGSRPDzb9gZeYwqS/jtuw/gJjs1lcfd17VuNPjgeblYwy07hrBQaXg8iapu9qy5XCe4dkse3zxiNYsO82VkiTieXNSs+ShktFYmwRpyRt0b/vkygJslxl4zl1CgysSRuYOMa8Bdb+I2v1wJBJOIACuK4Udzslx5v3k9W7I+LH++Mts9WO2JO2zwZ7SOu/RjOReHKkuhDf7CGvjxxVBxjGsASNoLYtR7YO81DPPl+rJMuaQqd/S3CZudcML+nPnMJh537LNmmHz96XMArMygwOwmw4QcUkvQb/zy7Vfh1l3D7Q+EGyQopbhcqGMin3RqZhiqDXPZDKmb4GdTh8lNgCv7LPezYK/NZzix+ybqmuNVBQZn9kXSnerHvMakKrXMdfe35pivxGt+uRKs/hW7hvGWA+P4+gNvwCHuhhpIqdg8kGwxr+dKdQyltRbdmm8DAMSjuG52k7lsShwFS24KLqZj5xlmXAPtZ0kwsEllUbvndnLTfEXvaXofQz6pBAYJthkIy/PfNZrBzpE0Li7VMGynPwYNU6rrTSdoricMpFUsVnUsVXU0jCbG80nkkgGexFoIEpvd2dRRxXusnX7cthx+sNfmA2WcOSsOkwgIErw8xqSjHcPplkJQVZYwlFY5JqFHJjN0A+vvqu0jZIng9RyLYNi/KRfAJFqrrQF3oWSdVeN6EoCVi75ceSIKaoRxHTSVDvBque2m0jEkFRYkIpiE7b+EYb5c76neysCa/PlxfKqIXELxSIt+3L7XYhPsGLarrvuYRNCc77UO9ne5bNdITOQTyCZVL5PQWycirgau4VqtRyVXsMU6blsOP4KCxJLDwNt7EpUAuYlnPm6QCE7fHs0mMFu0vIi5PtwfIkgsA/sncjg9XfJUTs6VWqutAXfXwnYaceUm6zlGj5gECTWunYu9xZPgmERsT8KWmyI8CZbJFYb5UqOnOeAMYcb18ali6MQ4BiY5uUHCbUnCUFuvTCJl9QebsmskJhwm4a2TSK6BIDGSTTiFb1FMIrVCuSkou4ltriLlJq01640FGl7OHeKYRBDGcgk3u6ms9zSzCRBBYlnYvymHhtn0zNWdLdUD0zTZhcEWoLjzJACrYKkXTMLq3RRmXAd7EoosOecSJwUWcN/7yuSmRl/kJl57Z6CU4thUoW1LiTdePYKkKjn9/oNSYOtGc03o9p1iMKXCaFK8ag/1mcglrbRYewdsNikaRnNNyE2A5Uuk7ZTmMKTt2RnLlZuYJ1equ59vHLnJSYHlmEQxQG5i13tY+rbVmsP1JHrNJER20zLAFo2XLxWd2dGhTMKRm2xPImZ2k/UcA9tDOkquBFHFdFEGXCahoG40YqXAAjHlJiVcbqo0DNT0ZtvZ1t3AQEqF2aQoN0xnF3ppqYZCzcCBNkEik1DwyMfehAl7F6tIBBLxehI1vTf+Uq/BroMT9njQ8XwCuaTieGaG3ShyLchNAPA/HNrWtjZmpXJTQpEgS8TTCTfOxEZ2P/DZTeV6q9zEmEEYkxjNWq05KKVYKPeeaa+/q3YNYP9EDmlNxjN24U5NN1GsGy0tOQBvy2wgXlOxBMckemVch7XlKNR0yBKJ7O3TSQos0IZJqOHZTW61de/TK4Oqro87pnXwxDge+yZyzmsQQloGzNT0tZEB1CmcIDFVxGBaRVKVXU2+Zrjdh9fIe7v7us349LuvjTxmpdlNhBC7yZ/7+YbJtDwke6ysx7i2W+/zfz+2joTVxIzlEqg0TFwu1GE0ac+ZtmASy4AiS7h5xyCePmtNQAsrpANcw3c5nkS9TZ3EcqFG1EksVXXkk0qgBs8Wh6j0Qh5OkIj0JCSPwcuDFdL1g0lM2O1MXp0pO95Cu8ymKCRV2de7qTefZa/B+lqduFx0qrGZ3FiqG5Ds6ySoOG+twmESK/g8sr6q80LNSlltl7Lun3PN+jbx99u7btgCWZJC+46xIMIyLIUnsUZx665hvDJVQKGmO7OO/R1gAXeHsORkN8WbJ8HQdyZRNULNt7QmI63JsXdgjnHdRm4KO5e5PlSTMrx+9zA0RcLjx6edx45PFbBlILmsWpSkInlbha9zJlGoGU4HArZJKNYMZ1e8VphEHKRXKDcBVkq417gOnwjo/d2K15PgptIxDGU0/I+v3xGaLMFqJ1iQWNPZTYSQYULIdwkhJ+3/h0KOu98+5iQh5H77sTQh5NuEkGOEkKOEkD/mjv8FQsgMIeSI/e+XVnKevcCtu4ZBqTUBLazaGuD9hc49Cevr7t987YzrMPMtk1BiF9IBaNuWw/pZuNw0v4zmfstFWlPwxqtG8PgxN0gcszObloOET26yjOv1tyfjAyTzXHJckHCHXa2f98bY/UqZhL9OIk79UFKVPNdFOWIedxiYYuEwibUcJAB8AsC/Ukr3AvhX+3sPCCHDAD4N4PUAbgPwaS6Y/Cml9ACAmwH8NCHkHdxTv04pvcn+9/kVnmfXcfOOQcgSwdNn550hIGMBTEKSrB5GnWQ38YGhJ9lNsgQ9ZDJd1I5oLOemGMZBPE8iPLvJlZt6HyQA4K0HxnFmtoxXZ8vQzSZOz5Ri+RFB8BvyNd3sScDvNfiGdWyqHbs+ijWdm62+juQmdWWeBGBtmPx1EnEYJz+EDPC2CY8LNjXzuD3Tu5dT6YCVB4l7AHzJ/vpLAN4bcMzbAXyXUjpPKV0A8F0Ad1NKK5TSxwGAUtoA8CyAbSs8n74hrSm4bkseT5+NZhKAdWG4nkR849r/dbegyhLMJoXZbA0USxFNyj71roN48MOHYv8et05iedlNc+UGFInELt5bKd56YBwA8NixaZyZKUM3qadAqxO0Gtfrk0lkuHTSCV8NgiU3WQvlaneB7QRMblpu7ybAYlMlX51EnOs0rSqeeRLFutHxaNfhjAZCgFMOk1jbFdcTlNJLAGD/Px5wzFYA57nvJ+3HHBBCBgG8GxYbYXgfIeQFQshDhJDtYSdACHmAEHKYEHJ4ZmYm7LCe4NZdwzhyfhGXFmu2Xh8+DavgZDd1Kjf1wri2bo4g87pQC/ckBtOas1DEQSwm0UZuGspobWdUdwvbh9PYO57FY8cu49iUNQpzuXJTUvU2LrTavq+fhZSBdcgF3Nb5bp2Agao9G2E9eRIsoK3E78to3hGmseUmTUaVuy7KEVP0wqDIEobSGsoNE5os9XR+OBAjSBBCvkcIeSng3z0xf0fQHe5sYQkhCoCvAvgspfSM/fC3AOyilN4A4Htw2UrrC1H6IKX0EKX00NjYWMxT6g4O7RpGw2ji8ePTkfMOUprsVGTGnSfB0CvjGgie4xDXgIuDOEGCdeT0z+cA+ldIx+OtB8bx1KvzeOa1BSgScWZCdwp/dtN6ZRKAVVAHuHKTa1zrjnSyVuok4oB1XI1zL4Yhm1RaGvzFMq5VGVWOSZRqnTMJwJW2hzJqzzdRbf9KlNL/QCm9LuDfwwAuE0I2A4D9/3TAS0wC4JnANgAXue8fBHCSUvrfud85Ryllg1z/CsAtnb2t/uDQLstamVyoBqa/MrDBQ0Dn2U29Mq6BViZR003UjWbsBn7tcHBLHleNZjAewT6iptPNlxs9T+/z460HxqGbFA89M4mrx7LLDtIJxTUoKaWoGeszuwlwq4gZi0woMjRFQrFuOAveenpvd107gT+459pYA6nCwI8wpZRGMnAeKU321Fcsx5MAgNGcPc62D/fHSrc2jwC43/76fgAPBxzzHQB3EUKGbMP6LvsxEEL+EMAAgN/gn8ACj433AHhlhefZE4xmE7jKHiw/EpHLz99AHXsSPdh9sh2UP0jEaS3QCV63YwiP/e93Ru6U2gWJ4Q4m0nUDt+wcQj5ppSkuV2oCrM+cvSfdpKC0N9JhPzCQUkEIPG2rc3adgJMCu46YRC6p4kNv3LWiHTg/wrTcMGE2aaxOBLtGMri4VMVcqY5mk6JUb02BjQPGJPqR1LHSq/aPAfwMIeQkgJ+xvwch5BAh5PMAQCmdB/AHAJ62/32GUjpPCNkG4HcBHATwrC/V9eN2WuzzAD4O4BdWeJ49w607rTbiQdXWDN65uh3KTT2YP8DOwV+f4HSA7ZNRDLimdlCG03y5/3KTIku4Y79lrR1YpmkNWC0YGJNgstN62m3zGM5oGM0mPNduzpZb1qMn0Q3wTf7itORguHP/GCgFfnhyFhX7+ug0BRZw02D70fxyRasBpXQOwNsCHj8M4Je4778I4Iu+YyYR7FeAUvo7AH5nJefWL9y6exhfP3w+Um7ib6A48oXWayYRIjd1m0nEgcMkfJ1gpws1LFV1J92vn3jrgTF86/mLuGaZ6a8Ay4e33pMzlW6dLqS/9par8b7XeRMPs3Yn2IpudSru5dCbtQh+hCm7j+LcN9dvHcBIRsPjx6edYWbZROf3G2N1vU5/BURbjhXj9buHQQiweTBcd091yCSsCWYEukmhyd1fWNQQ43qpgx1RtxAmN33lKWvS2ztv2NK3c2F41w1b0GwCt+9bfiIEnwLL/l+vctOe8ZzTyJIhl7DmXNcaa2OWRL/BjzBltSJxPAlJIrhj3xgePz6N//lOa544P9ArLkbXkdy04bF9OI2//+hPtey0eKQ69CQAV3LqxcKScJiEN6MozgjG7p9Lq9ykm0185clzuH3fGHaH9K/pJVRZwvtu2bai3TGr/6CUOu9tvcpNQcjadQJrZSpdv5Hhmhw6c+Fjbq7uPDCOhYqOH5+eBRB/HDCP0ZwIEusKN+8YilwAvEEi3p+cSUK9kJtCPQl2scdsBd4NsPfHM4l/OXoZ08U6PvyGnX07j27D9Vqajuy0XplEENgI00pjYwYJvlakEKMDLI/b945CIsA/vnAJwPLkJtaEspO6peVCyE19QKfGNeAuKL0xroOL6Tox4LqFIE/iy0+cxdbBFN5yIKg2c33AHWHavCKZRD6poljTUdPNdZXZ1C04TKJuhM6FD8NgWsPNO4bwzGtWF+nl1EnsGc/iG7/yRtyyM7BdXldx5Wxt1jDYaEdZIrEljITDJHpXJ+Fv8leoWpPw+rmY+eWmE5eLePLVefz8G3auazPUGWFqmA6TWI/jS8PAGtyV6xuTSWS57Cbm5XWy2N/J+V3LrZi+bfdwX+6RK+eqXcNIOw3F4n+gbPHsawpszNYC3QQLhmwh/dsnXoOmSPi5W0M7sawLsM/PKlC0jesraDHNJRU0qZWmvBGZhDN4qW6iULWqppUO7lWeJS8nBbafWNtnd4VgOYPXE2rvPIlEWApszCZl3QST4n71/3kGCXsY0ntv3tq3zq+9QpLzWhwmsU7bcgSBLWzTxdqKKpfXK5KqO8K0UNM7vm8Obs47s6qXk93UT4gg0QekltErpreexNphErtHM/jjn70eU4Uaqg0TDbOJX/zp3X09h14gGcQk1mGDvzCw6XQLFX1DMglCCDJ2i41CtfP7RpII3rJ/DI++eGnNXxciSPQBy+lf76TA9rGYbqmq930HTwjBvbft6Ovv7AeYr1PTr0wmwbeS2IieBGAFymLNsNvrd765+u13HFgXsqoIEn2AMy5R6cSTaD/6c7kIb8uhY9dI/+sSrkQ4xrVuOsV0yTW+Y+wEfG7/RmQSgDvCtFAznJTUTjCaTUR2algruHK2NmsYyWUwCadOohetwp3sJl8xXc3oa43ElQw3a6vp1ID0ghWuFnizdaMyCZbhVYgY1HUl4Mq5atcwGJNYK56EFsAkKKUoxBzBKNAeVz6TcK+TjRok2AjTQu3Kvm9EkOgD2E3UyWyChCJDlQmkHuRBBxXTVRomjCbtayHdlQzXk7BmdGiy1JPPcrXA5/ZvVLkpl1RQqOko1owr+r4RQaIPWG4KbK+yHhRZgkS8QWI1OsBeyUg4xXRN1HTzimrJAYggAVhN/qaWagCu7PvmyhXS1hDcIBF/J/mBW7Zj38Ty5xm0gypLHrnJnSVx5V7s/YTblsOquL6SCukAq3sA0+Q3YhdYwJKbKh10gF2vEEGiD1hOCuz12wZw/baBXp2SNVs6gElcyRd7P8G3QK8bVx6TAFzjdqN6EnyGV7+LUPuJK+/KXYNQZQmqTHpiQi8Xmo9JLFU662QpEA1NlkCI7UnozSuqRoKBLZJXUuPCTsDPpr6S5aYr78pdo0iqckdMotfQ7BYYDI4nIeSmroAQ4owwrenmFbmQsjTYtLYxNxa8L3Ml3zdrZ9W6wpHWZKhrSHJo9SSE3NRtsBGmdaN5RcpNLA12o8pNfJAYSF+5982KrlxCyDAh5LuEkJP2/4HNzf//9u4+RqrqjOP497fzsqu8CMirQKvGTX2pEe1W7UsaqgbRGrGpVo21q4EYExNt06bBmoZUY6JJU23TxgQRi01LNWgq7R81iJr6T4lLNdbWGoitsHaBlYWVsMCy7NM/7pndu+tcdpeZcZh7n0+ymblnztw5h3OZZ84599wrqT3k2SapPZb+mqT3JL0V/maH9GZJz0raLmmLpDMrKefJ4ML50zh3bu0moieqdHvUktINh07kLlmuvJZCjiMD6e1JlC7NcUoxfQFwPEb2JNL7/6bS1l0JbDazVmBz2B5B0gxgFXAZcCmwalQwuc3MFoW/PSFtObDPzM4BHgMerbCcdbemvY17vn5OvYsxpJjPjZi47j10lEnF3IQud+yOrzk/3JNIZZBoKQWJ9H5BHk9pTqJJw/e8TqNKvxGWAevC83XADWXyXA1sMrMeM9sHbAKWTmC/G4ArJaVnJdJJoJjTJ4ab0jz5Vg8theE5iTQON5V+SWd1uKkUJKeeUkjVQsnRKj1y55hZF0B4LHe/yfnAzth2Z0greToMNf0kFgiG3mNmA0AvcHq5Aki6S1KHpI7u7u7KapMh5SaufT6iupoLuWgx3UA6h5umnVpAItPrJCDdk9YwjnUSkl4G5pZ56YFxfka5EFsaDL/NzD6UNAV4HrgdeGaM94xMNFsNrAZoa2srm8d9UrnFdGk/2D9tLfmmoVNg09iTuPmLn6F1zpRUBsDxKPWk0n7a+Ji1M7Orkl6TtFvSPDPrkjQP2FMmWyewOLa9AHgt7PvD8HhA0u+J5iyeCe9ZCHRKygOnAT3jqZAbn2K+iYNHBoa2ew8d5YxpLXUsUfo0F3L0Hjqa2onrWVOaufqCcr8fs2FyRnoSlf682QiUzlZqB14sk+clYImk6WHCegnwkqS8pJkAkgrAdcA7ZfZ7I/CKmXkvoYoKuaahS1hDuCtdyg/2T1tLvim6LMfAYKouE+4ipVuYpn2YttJ+0iPAc5KWAzuAmwAktQF3m9kKM+uR9BDwRnjPgyFtElGwKAA54GXgyZDnKeC3krYT9SBuqbCcbpT4nISZ0X3gCDOnnPw3QGkkLYUch44eo39g8KS/RaWbOCm6flXaf1xVFCTMbC9wZZn0DmBFbHstsHZUnoPAFxL2e5gQcFxtFHNNQ+sk9vcd5cjAIHOn+nBTNbUUmoYWKabxshwOHlx2Aa2zT571T7WQ7hkXl6gQOwW2K1zueN5pHiSqqSXMSUBtbkPr6m/ZovljZ2pw/vMmo+LDTbs+PgTAXA8SVdWcb2IwzKR5T8I1Kj9yMyp+CuxwT2LiN3N3yeJnNKXp1qUuWzxIZFT8fhJd+w+TaxKzfOK6quJBws9uco3Kj9yMKuaiIGFmdPUeZs6UZnIpvrRAPcQX0HlPwjUqDxIZVcw1YQbHBo1dHx/y+YgaGDHclMLFdC4bPEhkVOneFv3HBunqPezzETUQ70n4cJNrVH7kZlTpVqr9A4Ps6j3sPYka8IlrlwYeJDKq1JPYe7Cfvv5jvkaiBnzi2qWBH7kZVcxFk9Q79vYBvkaiFuJrI7wn4RqVB4mMKoaexAd7DwK+2roW4qusfTGda1R+5GZUIcxJfNBT6kn4xHW1xQODX5bDNSoPEhlVmrje2dOHBLN9IV3V+ZyESwM/cjOqMDTc1Mesyc1DPQtXPfF5iDTemc5lgx+5GdUcgsKOnj6fj6iR0nBTc76J4du3O9dYPEhkVKkncWRg0BfS1UhpHsJXW7tG5kEio4qx4SU//bU2mmM9CecalR+9GRWfg/DhptqIhpm8J+EamweJjCrmh8fIvSdRG5Jozjf5GgnX0Co6eiXNkLRJ0rbwOD0hX3vIs01Se0ibIumt2N9Hkh4Pr90hqTv22opy+3Unrpgb/nXrcxK101LI+RoJ19Aq/YmzEthsZq3A5rA9gqQZwCrgMuBSYJWk6WZ2wMwWlf6AD4AXYm99Nvb6mgrL6UYpxHoSPtxUO96TcI2u0qN3GbAuPF8H3FAmz9XAJjPrMbN9wCZgaTyDpFZgNvB6heVx4xSfuJ491RfS1Yr3JFyjqzRIzDGzLoDwOLtMnvnAzth2Z0iLu5Wo52CxtG9JelvSBkkLkwog6S5JHZI6uru7T6wWGVQ6BXbm5KJ/idXQ5OY8k5r939c1rvxYGSS9DMwt89ID4/yMcquIbNT2LcDtse0/AevN7Iiku4l6KVeU27mZrQZWA7S1tY3er0tQ6kn4pHVtPfzNC5nsQcI1sDGDhJldlfSapN2S5plZl6R5wJ4y2TqBxbHtBcBrsX1cBOTNbGvsM/fG8j8JPDpWOd3ElE6BnTvVJ61radHCafUugnMVqXS4aSPQHp63Ay+WyfMSsETS9HD205KQVnIrsD7+hhBwSq4H3q2wnG6UXJPIN8knrZ1zxzVmT2IMjwDPSVoO7ABuApDUBtxtZivMrEfSQ8Ab4T0PmllPbB/fBq4dtd97JV0PDAA9wB0VltOVcf+153HZWTPqXQzn3ElMI+eKG1tbW5t1dHTUuxjOOddQJG01s7Zyr/kJ3M455xJ5kHDOOZfIg4RzzrlEHiScc84l8iDhnHMukQcJ55xziTxIOOecS+RBwjnnXKJULaaT1E10X4oTMRP4qIrFaRRZrHcW6wzZrHcW6wwTr/dnzWxWuRdSFSQqIakjacVhmmWx3lmsM2Sz3lmsM1S33j7c5JxzLpEHCeecc4k8SAxbXe8C1EkW653FOkM2653FOkMV6+1zEs455xJ5T8I551wiDxLOOecSeZAAJC2V9J6k7ZJW1rs8tSBpoaRXJb0r6Z+S7gvpMyRtkrQtPE6vd1mrTVJO0puS/hy2z5K0JdT5WUnFepex2iRNk7RB0r9Dm38pI239/XB8vyNpvaSWtLW3pLWS9kh6J5ZWtm0V+WX4bntb0iUT/bzMBwlJOeDXwDXA+cCtks6vb6lqYgD4gZmdB1wO3BPquRLYbGatwOawnTb3MfI+6Y8Cj4U67wOW16VUtfUL4C9mdi5wEVH9U93WkuYD9wJtZvZ5IAfcQvra+zfA0lFpSW17DdAa/u4Cnpjoh2U+SACXAtvN7H0z6wf+ACyrc5mqzsy6zOzv4fkBoi+N+UR1XReyrQNuqE8Ja0PSAuAbwJqwLeC+Cr96AAACYElEQVQKYEPIksY6TwW+BjwFYGb9ZraflLd1kAdOkZQHTgW6SFl7m9lfgZ5RyUltuwx4xiJ/A6ZJmjeRz/MgEX1R7oxtd4a01JJ0JnAxsAWYY2ZdEAUSYHb9SlYTjwM/AgbD9unAfjMbCNtpbO+zgW7g6TDMtkbSJFLe1mb2IfAzYAdRcOgFtpL+9obktq34+82DBKhMWmrPC5Y0GXge+J6ZfVzv8tSSpOuAPWa2NZ5cJmva2jsPXAI8YWYXAwdJ2dBSOWEcfhlwFnAGMIlouGW0tLX38VR8vHuQiCLrwtj2AuB/dSpLTUkqEAWI35nZCyF5d6n7GR731Kt8NfAV4HpJ/yUaRryCqGcxLQxHQDrbuxPoNLMtYXsDUdBIc1sDXAX8x8y6zewo8ALwZdLf3pDcthV/v3mQgDeA1nAGRJFoomtjnctUdWEs/ingXTP7eeyljUB7eN4OvPhpl61WzOx+M1tgZmcStesrZnYb8CpwY8iWqjoDmNkuYKekz4WkK4F/keK2DnYAl0s6NRzvpXqnur2DpLbdCHw3nOV0OdBbGpYaL19xDUi6lugXZg5Ya2YP17lIVSfpq8DrwD8YHp//MdG8xHPAZ4j+k91kZqMnxRqepMXAD83sOklnE/UsZgBvAt8xsyP1LF+1SVpENFlfBN4H7iT6UZjqtpb0U+BmorP53gRWEI3Bp6a9Ja0HFhNdDnw3sAr4I2XaNgTLXxGdDdUH3GlmHRP6PA8Szjnnkvhwk3POuUQeJJxzziXyIOGccy6RBwnnnHOJPEg455xL5EHCOedcIg8SzjnnEv0fAtX13ZWbw8AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "lambda_F = 8e-4\n",
    "sigma_F = 3e-2\n",
    "sigma = 2e-2\n",
    "\n",
    "# Create the factor time series\n",
    "F = np.random.normal(lambda_F, sigma_F, 100)\n",
    "\n",
    "returns = {}\n",
    "\n",
    "# First group\n",
    "beta1=0.8\n",
    "for idx in range(10):\n",
    "    returns[f'stock1_{idx}'] = beta1*F + np.random.normal(0,sigma,100)\n",
    "\n",
    "# Second group\n",
    "beta2=0.3\n",
    "for idx in range(5):\n",
    "    returns[f'stock2_{idx}'] = beta2*F + np.random.normal(0,sigma,100)\n",
    "    \n",
    "# Third group\n",
    "beta3=0\n",
    "for idx in range(5):\n",
    "    returns[f'stock3_{idx}'] = beta3*F + np.random.normal(0,sigma,100)\n",
    "    \n",
    "plt.plot(returns['stock1_2'])\n",
    "\n",
    "# Dump into array\n",
    "returns_arr = np.array([r for r in returns.values()])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1a170ac588>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAD4CAYAAADB9HwiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de5hdVX3/8ffnzGSSTO4QgXAT0IhYqgEpRa2gcikiBWmlhao/rFZaK4q0/iz+6INWnz6lqLVa++jDTfAGKqIiRQVRpPYR5I7BhIsUIRASLobcL5P5/v44m/Qwmctea++Z2TP5vPLsZ86c2d+zVs6cWWedtdd3LUUEZmbWPK3xroCZmQ3ODbSZWUO5gTYzayg30GZmDeUG2sysobrHsrCXvnv35Ckjbzl5YXI5Tz6zJjkGoKsr/f1q/txZyTHTe6YnxwCs2bA2Ky7H+o2bk2N6pqS/nM444pPJMQBPbXg4OWZz/8assnLsNG2v5Jj9Nh6aVdbq2b9Nj+l7NKusaV1zkmN22rRvVlk982YoK7BD2TZn6YWPVy5rNLgHbWbWUGPagzYzG0ua4F3QStWXdKyk+yQ9KOnsuiplZlaHri6VOpoquwctqQv4D+BoYBlwq6SrI+JXdVXOzKyK1g7cgz4UeDAiHoqIzcAVwIn1VMvMrDq1VOpoqioN9B5A56XgZcV9zyPpdEm3Sbpt1dL1FYozM0vTapU7mqpK1QZ729luSktEXBARh0TEIXNf2luhODOzNGqVO5qqyiyOZUDnZM89gcerVcfMrD4tNXf4oowqDfStwEJJ+wKPAacAf15LrczMatA1wScSZ1c/IvoknQH8EOgCLomIe2urmZlZRU0eviij0vtLRFwLXFv2/Jy07Su/+UByzHHHvzA5BvLSm1ev25Acs3Z9XsrxzN5pyTGtzCvUOf+vnI+T67akpykDrNo0dqNpU7tmJsc8szE9lfrFG1+ZHAPQPTtv6YAc87bulxzT2rJ1FGpSsuwGz9AoY4J/ADAzG9oEH4J2A21mk1eTp9CVUTXV+xJJKyUtrqtCZmZ16epWqaOpqr6/XAocW0M9zMxqtyPPgyYibpK0Tz1VMTOr10SfBz3q7x2dqd53XLt8tIszM9tmovegR71qnaneBx+3YLSLMzPbZkdei8PMrNHqXM1upPXvJe0t6SeS7pR0j6Tjqtbf0+zMbNLq6qrncUquf/8PwDci4vOSXkY7iW+fKuVWnWZ3OfBzYH9JyyS9q8rjmZnVqdVSqaOEMuvfBzC7uD2HGhaPqzqL49SU83N2285J2772mt8kxwCc8e7XJcc8seqx5JgZU2ckx+Ra/vTTWXE93ekvjZ3mzB75pAFayuvi5KRf90V6Kj9AV6snOWZ6d/pzEVv6k2MAlNHPmt2dvus4QGTMiujvG7+R1LIXACWdDpzecdcFEXFBx/eDrX//+wMe5qPAdZLeB8wAjkqs7nY8xGFmk1bZtTiKxviCYU4ps/79qcClEfEpSa8CvizpwIjIe+elwhCHpL2KAfElku6VdGbuY5mZjYaWVOooocz69+8CvgEQET8HpgHzK9W/Qmwf8HcRcQBwGPDeYmDczKwRurtbpY4Stq1/L6mH9vr3Vw845xHgSABJB9BuoJ+sVP/cwIhYDiwvbq+RtIT2OI139TazRmjVlIUy1Pr3kj4G3BYRVwN/B1wo6Szawx/viIjttgFMUcsYdJHufRBwSx2PZ2ZWhzrXgx5s/fuIOLfj9q+A19RWIDUkqkiaCXwL+EBErB7k59tSvZfc8FTV4szMSqtxmt24qDoPegrtxvmrEXHVYOd0pnofcGSl8XIzsyQ1XiQcF9lDHJIEXAwsiYh/ra9KZmb1aDV5oY0SqoxBvwZ4O/BLSXcV9/2/YpzGzGzcdXftoA10RPyMwSdvD6kr48nK2cg1JyMQ4HMX3pgck5PpuH7DZvbeZc/kuC1btyTH9EzJy9TL6Xls3pL+u/rqrf/MyQd/IDluald6Nubsrl2TYwCe2ZS+Aey8qbunF9Sb9+copedBPLbuzqyyFvb+YVbceGny+HIZziQcBzmN82SV0ziblbUjD3GYmTWaGnwBsIwqFwmnATcBU4vHuTIiPlJXxczMqtqRhzg2AW+IiLXFdLufSfp+RNxcU93MzCrZYRvoIoVxbfHtlOKolNZoZlan7rpW7B8nVRNVuoopdiuB6yNiu1TvzkzCX/2o0rohZmZJJnqiSqUGOiK2RsQi2kvvHSrpwEHO2ZZJ+LKjXlClODOzJBM91buWWRwRsUrSjcCxwOI6HtPMrKqJPs2uyoL9L5A0t7g9nfb2LkvrqpiZWVUTfYijSg96AXBZsdtti/ZuttfUUy0zs+pKLsbfWFVmcdxDew3o0ubPnZVczup1G5JjcjZyhbHboPbNJ+U97dOmTkmOyUmVB5jWk15WX0Yq/6ata0c+aRDL1z6QHLNg5sKssmZ0z0uO6evPeN4z13bfEumbMb+w99VZZXVt6EsPyt6Rr7q6FuwfL84kNLNJq8kXAMtwA21mk1ZLE3sedOUGuhiDvg14LCKOr14lM7N6TPRZHHX0oM8ElgCza3gsM7PadE3wHnTVTMI9gTcBF9VTHTOz+nR39ZQ6mqpq///fgA8xzHXazlTvO/7z8YrFmZmV12q1Sh1NVSVR5XhgZUTcPtx5naneB78pY5cJM7NMXeoqdTRV1T0JT5B0HDANmC3pKxHxtnqqZmZWTavV3Ma3jOwedER8OCL2jIh9gFOAH7txNrMmaalV6mgqz4M2s0lroveg61rN7kbgxpHOm94zPfmx167fmBwzY2r6js8A82ftkhyTk7b9nW//OjmmXdaLkmN2f8HOWWX1Zewg3rd1a3LM9O45yTEA83v3So5Z1/fbrLJ2631JcszmrelLFJC5aM90pS/ju6r/wayyujKe954149dDndJKX7KgSdyDNrNJyz1oM7OG2qFTvSU9DKwBtgJ9EXFIHZUyM6tDl3vQvD4inqrhcczMatXkGRpleIjDzCatJqdxl1H17SWA6yTdLun0wU7oTPX+xfceqVicmVl5Ez3Vu2oP+jUR8bikXYDrJS2NiJs6T4iIC4ALAP75p8flbRlhZpahyWncZVR664iIx4uvK4FvA4fWUSkzszq01FXqKEPSsZLuk/SgpLOHOe8tkkJS5UkTVRZLmiFp1nO3gWOAxVUrZGZWl1arq9QxkmJjkv8A3gi8DDhV0ssGOW8W8H7gllrqXyF2V+Bnku4GfgH8Z0T8oI5KmZnVoca1OA4FHoyIhyJiM3AFcOIg530cOB9IT4EeRJVdvR8CXpESs2ZD+g7OM3unJcfk2pKR3pyz03ZOyjbkpYj/xf95ZVZZc3vTU8Q39aW/JudP2y85BmDGlPSdtqe00pcaAFi7JX0W6Ya+1ckx/VNz+0vpl3a2Rt5u7zlxMaU3q6w6TKlvFscewKMd3y8Dfr/zBEkHAXtFxDWSPlhHoZ5mZ2aTVtl50MUstM6ZaBcUExy2nTJI2LZ3Rkkt4NPAO9JrOTQ30GY2aZW9ANg522wIy4DOlaL2BDq3iJoFHAjcqPaiV7sBV0s6ISJuS6lzp6p7Es6VdKWkpZKWSHpVlcczM6tTjbM4bgUWStpXUg/tNfCvfu6HEfFsRMyPiH2KNfJvBio1zlC9B/0Z4AcR8Zai0uM32GRmNoBqmgcdEX2SzgB+CHQBl0TEvZI+BtwWEVcP/wh5shtoSbOBwynGXIorm3lXHszMRkGdq9lFxLXAtQPuO3eIc19XR5lVhjj2A54EvijpTkkXFfOhn6cz1fvOa5dXKM7MLE23ekodTVWlge4GDgY+HxEHAeuA7bJrOnf1Pui4BRWKMzNLU2cm4Xio0kAvA5ZFxHMZM1fSbrDNzBpB6ip1NFWVXb2fAB6VtH9x15HAr2qplZlZDVp0lTqaquosjvcBXy1mcDwE/EX1KpmZ1WOHXrA/Iu4CRnWbq1Yrfafj5U8/nVVWz5T0d9L1G9MnruTutJ2Ttv3FL92eVdYpp7w0OeaZZ9clxzy18aHkGICHVqVPL120y5uyylq88sfJMQftdlxyTGt1+q7oAKumPz7ySQPMaO2aVVYPs5Jj+sdutYbtdLeaewGwDGcSmtmk1eTx5TLcQJvZpNXk8eUyqqwHvb+kuzqO1ZI+UGflzMyqmOjT7KosN3ofsAi2LWb9GO1dVczMGqHJjW8ZdQ1xHAn8OiJ+U9PjmZlVNtEb6LrmoJwCXD7YD5zqbWbjpUs9pY6mqtxAF3OgTwC+OdjPneptZuOlxi2vxkUdQxxvBO6IiBU1PJaZWW0m+hBHHQ30qQwxvGFmNp526AZaUi9wNPBX9VTHzKw+muDzoKumeq8HSuct56RFr163ITmmpzvvv9VqpY9FTetJ39W7L2P3cMjbaTsnZRvgiiuWJsec9rZFyTFb+tN/vwB7zf6d5JhH1tyZVdb0nu2WOR85pntOcoy29CfHAExT+utiY2Quh9CanRzTvSFzjLeGFHGnepuZNZRqm6g2PtxAm9kklr7YWpNU3dX7LEn3Slos6XJJ47hulZnZ84lWqaOpqqzFsQfwfuCQiDiQ9k63p9RVMTOzqlTyX1NVHeLoBqZL2gL0AukL05qZjZrm9o7LqLLl1WPAJ4FHgOXAsxFx3cDzOlO9f3ndyvyampklmuir2VUZ4pgHnAjsC+wOzJD0toHndaZ6/+4xu+TX1MwsmUoezVSl/38U8D8R8WREbAGuAl5dT7XMzKqb6BcJq4xBPwIcVmQTbqC95Gj6RnFmZqOkyRcAy6iyYP8tkq4E7gD6gDuBC+qqmJlZdc3tHZdRNdX7I8BHyp7fMyW9uJbS3wF3mpOejgqweUt6KnpfV/oLoG9r3u7Nm/o2Jsfk7LQNeWnbl33lruSYN/5zcggAT214ODlmj5kHZpX17Kb0hRrvX/VfyTE7z/2L5BiAzaRffH903R1ZZR045Y+SY2LK+DWSO2wP2sys6Sb6YklVMwnPLLII7/WGsWbWNBP9ImGVaXYHAu8GDgVeARwvaWFdFTMzq2qiZxJWees4ALg5ItZHRB/wU+CkeqplZlaHVsmjmarUbDFwuKSdi6l2xwF71VMtM7PqdtghjohYAvwLcD3wA+Bu2tPtnqcz1fvuHzyRXVEzs1RSV6mjqSq9dUTExRFxcEQcDjwDPDDIOdtSvV9x7G5VijMzS7Ijj0EjaZfi697AH+PNY82sQeoc4pB0rKT7JD0o6exBfj5V0teLn98iaZ+q9a86D/pbknYGtgDvjYjfVq2QmVl96ukdqz0O8h+0N8leBtwq6eqI+FXHae8CfhsRL5Z0Cu0h4D+rUm7VTMLXVok3MxtNNV4APBR4MCIeApB0Be3VPDsb6BOBjxa3rwQ+J0kREbmFjmkm4RlHfDI5Zt2W9E557vqum7em7zC9aeva5JicHZ8B5k/bLznmqY0PZZWVs9t2Ttr2KR8+Mz0IuOkzP0mOCfJS7Oftkj45qVe7Jsd0/TavfrNIX8Z375mHZJWl9eltzcapq7PKmsr0rLhOZceXJZ0OnN5x1wUR0bm20B7Aox3fLwN+f8DDbDsnIvokPQvsDDyVWO1tnOptZpNXlGugi8Z4uMXeBnugge9WZc5JMmL/X9IlklZKWtxx306Srpf0QPF1XpVKmJmNBkWUOkpYxvPzPPZk+y3+tp0jqRuYQ3t2W7YyAzSXAscOuO9s4IaIWAjcUHxvZtYsUfIY2a3AQkn7SuqhvUH21QPOuRo4rbj9FuDHVcafoUQDHRE3sf27wInAZcXty4A3V6mEmdmoqKmBLpazOAP4IbAE+EZE3CvpY5JOKE67GNhZ0oPA31JDxzV3DHrXiFgOEBHLn5sPbWbWKNU6sAMeKq4Frh1w37kdtzcCJ9dWIGOwSkhnqvfXLvnOaBdnZraNotzRVLk96BWSFhS95wUw9JYOnVdHH1l7c4OfCjObdPondpOT24PuHAw/DfhuPdUxM6tRfRcJx0WZaXaXAz8H9pe0TNK7gPOAoyU9QDv18bzRraaZ2Y5nxCGOiDh1iB8dmVpYzkafqzYNnGo4sqldM5Nj2nEzkmOWr91uAb8Rze/NWzZ7xpT06eYPrbotq6y9Zv9OckzO7zcnIxDg8DNfnxxz42euyyorR06K8dp567PK2kJ6Nms3vVllRSt9bYv+zAzOOpSc49xYziQ0s8lrYrfPbqDNbBKb7BcJh0j1PrnYybtfUt6qK2Zmo2yiT7PLTfVeTHuB/pvqrpCZWW0m+CyOMhcJbxq4M0CxHyFSc7eKMTNrcuNbxphmEl71pRtGuzgzs/8VUe5oqFG/SNiZSXjHk1c095kws0mnyePLZXgWh5lNXg3uHZfhBtrMJq+J3T7npXpLOknSMuBVwH9K+uFoV9TMLNVEn2anigv+J7l5xaXJha3fkr5jTFerJzkGYKep6SnYOZvGrutL3wgX4MVzXp0c0x95abaPrLkzOWanaenP305TX5gcAxD0J8e87sxjssr62Wf/KyMq/e9qVt/uGeXAhu6811OO3o1zk2NaG/Neg90LZleeJtb/yKpSv4jW3nMbOSXNQxxmNnk1uHdchhtoM5u8dtBU709IWirpHknflpT+ucfMbJTF1ih1NFVuqvf1wIER8XLgfuDDNdfLzKyy6I9SR1Nl7eodEdcVu9wC3AzsOQp1MzOrZoJnEtaR6v1O4PtD/bAz1fs7X76xhuLMzMqZ6D3oShcJJZ0D9AFfHeqczlTvnGl2ZmbZGtz4lpHdQEs6DTgeODLGcjK1mVlJTb4AWEZWAy3pWODvgSMiIm8jNTOzUTbR+465u3p/DpgFXC/pLklfGOV6mpml649yR0M1PtV7a//m5HJyU71z0rZndKfvtD2r5wXJMbkWr/xxVtz0nvQdznMcvMtJY1IOQLemZ8X9wftfmxxz82fTd1OfuWZWcgzAptnpae9r+h/NKmuuXpwc070q/W8Y6kn13nj746XanGmv3N2p3mZmY2qCD3G4gTazSavJU+jKyE31/niR5n2XpOsk5S3DZWY2inbUVO9PRMTLI2IRcA1wbt0VMzOrbIJfJMzd1Xt1x7czmPCL+pnZpNTgxreM7FRvSf8k6VHgrQzTg3aqt5mNl4godTRVdgMdEedExF6007zPGOa8CyLikIg45M1vf11ucWZm6fpLHg1Vx2JJXwP+pIbHMTOrVWztL3VUJWknSddLeqD4OmSChKTZkh6T9LmRHjergZa0sOPbE4ClOY9jZjaaxnA1u7OBGyJiIXBD8f1QPg78tMyDjniRsEj1fh0wv9jJ+yPAcZL2p/3h4DfAX5cpzMxsTI3dRcITabeTAJcBN9Jer+h5JL0S2BX4AXDISA9aZhbHqYPcffFIcYPJ2fX5mY3pKanTu2cnxwDMm5o+nbsvIxV989YNyTEAG/pWj3zSAAftdlxWWdO75yTH3L8qfffrXu2aHAOgjA9/m8jb/Tonbfuw94/4t7edxR/L+yA6ZdOU5JjunmlZZSmjweubl14/qCeLbgwvAO4aEcuLMpdL2mXgCZJawKeAtwNHlnlQZxKa2eRV8g1F0unA6R13XVCsZd95zo+A3QYJP6dkbf4GuDYiHpXKLf3hBtrMJq2yFwA7NxYZ5pyjhvqZpBWSFhS95wXAykFOexXwWkl/A8wEeiStjYghx6uzUr07fvZBSSFp/kiPY2Y21vq3bC111OBq4LTi9mnAdweeEBFvjYi9I2If4IPAl4ZrnCE/1RtJewFHA4+UeAwzs7HX31/uqO484GhJD9BuF88DkHSIpItyHzQr1bvwaeBDDPJOYWbWBGO1EFJEPM0gF/4i4jbgLwe5/1Land9h5c6DPgF4LCLuLnHutlTvr196TU5xZmZZor+/1NFUyRcJJfXSvmp5TJnzOwff73/2huYmvZvZpFNHluB4ypnF8SJgX+DuYqrInsAdkg6NiCfqrJyZWSUN7h2XkdxAR8QvgW2TsCU9DBwSEU/VWC8zs8pqmqExbnJ39TYza7yJPgY9prt6961Yk1xYa2P6O2BsyXzCezNGfHKev5JZRAP1T02/ptvakNeDUMZz2Dc3fTf1rsz6rZ23PjlmWl/eEgBZdcyYPXDguS9NLwe4698fSo7pZ0tWWX2R/ry3NDWrrJnT5lbeaXvFF/671C9i179+jXf1NjMbSxN901g30GY2aU30WRy5u3p/tFhw+q7iyFsyzcxsFMWWraWOpspO9QY+HRGLiuPaeqtlZlbdRL9IWCXV28ys0Sb9EMcwzpB0TzEEMtz+W9tSvS/88hcrFGdmlmbS96CH8Hna+2pF8fVTwDsHO7Ez1Ttnmp2ZWbYxWixptGQ10BGx4rnbki4EvAqSmTVOk3vHZWQ10M/tHFB8exKw3WL+ZmbjrX9L33hXoZLcXb1fJ2kR7SGOh4G/GsU6mpnlmexDHHXu6r16dvquyt2zpyfH5Oz4DCClfxzaEmuSY6brBckxbekvtlXTH88qaZp2To7ZPOg2bMObxXabH5eyhbXJMdGdN991yuyZ6TEZO23npGwDLHrffskx95z/QFZZq6em/47nr3phVlksyAvrtEMOcZiZTQTR5wbazKyRJv086KF29Zb0Pkn3SbpX0vmjV0UzszzR11/qaKoyPehLgc8BX3ruDkmvB04EXh4RmyTlDSSamY2i/k2TfBbHEKne7wHOi4hNxTnpVw7MzEZZk3vHZeSmer8EeK2kWyT9VNLvDXViZ6r3ZRd/LbM4M7N0sTVKHU2Ve5GwG5gHHAb8HvANSfvFINuzdKZ6P7PhkeY+E2Y26Uz0HnRuA70MuKpokH+h9gTi+cCTtdXMzKyiiT6LI7eB/g7wBuBGSS8BegDv6m1mjdK/cZJfJBwi1fsS4JJi6t1m4LTBhjfMzMbTpO9BD5HqDfC21MJW9z2aGpJldvdeWXGPrbszOeaFva9OjlnV/2ByDMDW2JwcM6O1a1ZZG+Pp5JhH192RHLP3zEOSYwC66c2IyutDrOlPf91290xLjukl73eVk7b98g8tzCrr7k//Ojlm3dxns8qaQ94u7J121DFoM7PGm/QNtKRLgOOBlRFxYHHf14H9i1PmAqsiYtGo1dLMLMOkH+JgkEzCiPiz525L+hSQ9xnGzGwUTfoe9HCbxkoS8Ke0Z3SYmTXKpE/1HsFrgRURkbe4rJnZKJroPegqu3oDnApcPtwJnaneX/vidysWZ2ZW3o6a6o2kbuCPgVcOd15nqvfDa/67uc+EmU06O3IP+ihgaUQsq6syZmZ1iq39pY6qJO0k6XpJDxRf5w1x3vnFGvpLJH22uI43pDIL9l8O/BzYX9IySe8qfnQKIwxvmJmNpzFcsP9s4IaIWAjcUHz/PJJeDbwGeDlwIO2F5o4Y7kGzMwkj4h0jVtnMbByN4SyOE2kviQFwGXAj8PcDzglgGu21iwRMAVYM96Bjmkk4rWtOcsy8rek7FsfwnxqGtLD3D5NjujakvwC6evNS0XNSvXuYlVVWTys9zfbAKX+UHKP1eZclopX3O84xbVr6DufqT/9/bWqtTo6BvJ22c1K2AV5x1ouSY3J3K69D2eELSacDp3fcdUFx/aysXSNiOUBELB9sl6mI+LmknwDLaTfQn4uIJcM9qFO9zWzS6o9yDXTnZIahSPoRsNsgPzqnTBmSXgwcAOxZ3HW9pMMj4qahYnJTvRcBX6DdXe8D/iYiflGmkmZmY6W/xkU2I+KooX4maYWkBUXveQEw2Meak4CbI2JtEfN92pueDNlAl5nFcSlw7ID7zgf+sVh/49ziezOzRtka/aWOGlwNnFbcPg0YLOnjEeAISd2SptC+QDjsEMeIDXTR/X5m4N2wbS3AOcDjIz2OmdlY6+vfWuqowXnA0ZIeAI4uvkfSIZIuKs65Evg18EvgbuDuiPjecA+aOwb9AeCHkj5Ju5EfclHkzsH38//9XN7+zpMzizQzS1N2DLqqiHgaOHKQ+28D/rK4vRX4q5THzW2g3wOcFRHfkvSnwMW0E1e20zn4/sT6xc4kNLMxM1YN9GjJzSQ8DbiquP1N4NB6qmNmVp/+iFJHU+U20I/zvxkwbwC8mp2ZNU5/9Jc6mip309h3A58pFkzayPMneJuZNUJNMzTGjcZyM+7Nv12XXFhrS/oT3D+l6iqq5SlnqcLM5zwy/l/907qyymptSL+ynVO/TV152XP9pNdv5qqZWWXl6Js3JT2GTVllTftt+vO+bu76rLKmam5yzKL3pWcDAyy98PHK6aLfnPOeUn9sJz/7+bFLTU3gTEIzm7SaPHxRhhtoM5u0mnwBsIwyy41eImmlpMUd971C0s8l/VLS9ySlr6xjZjbKJvpFwtxU74uAsyPid4FvA/+35nqZmVU2hqneoyI31Xt//neBj+uBP6m5XmZmle0IPejBLAZOKG6fDAy5wHHnprEXXXpJZnFmZunGcC2OUZF7kfCdwGclnUt7FachV5LvTPXOmWZnZparyb3jMrIa6IhYChwDIOklwJvqrJSZWR0m+iyOrAZa0i4RsVJSC/gH2ov3m5k1SpMvAJaRm+o9U9J7i1OuAr44ajU0M8s06Yc4htrVG/hMzXUxM6tVky8AlhIRjTiA05saM1nLanr9/Fz4udjRj7FbVWhkOSvijVXMZC2r6fUby7KaXr+xLKvp9dthNKmBNjOzDm6gzcwaqkkN9AUNjpmsZTW9fmNZVtPrN5ZlNb1+O4wxXbDfzMzKa1IP2szMOriBNjNrqHFvoCUdK+k+SQ9KOrtkzHabCJSI2UvSTyQtkXSvpDNLxEyT9AtJdxcx/1i2vCK+S9Kdkq4pef7DxSYId0m6LaGcuZKulLS0+P+9aoTz9y/KeO5YLekDJco5q3geFku6XNK0kvU7s4i5d6hyhtgYYidJ10t6oPg6r2TcyUVZ/ZIOKRnzieL5u0fSt6XtN98bIu7jRcxdkq6TtPtIMR0/+6CkkDS/RDkflfRYx+/suDL1K+5/X/E3dq+k80uU9fWOch6WdFfJ52KRpJufe/1KOrREjDf/GM54TsIGuoBfA/sBPcDdwMtKxB0OHAwsTihrAXBwcXsWcP9IZQECZha3pwC3AIcllPm3wNeAa0qe/zAwP+N5vAz4y+J2DzA38e0JWs8AAAWVSURBVHfwBPDCEc7bA/gfYHrx/TeAd5R4/ANpL0/bSztz9UfAwjK/U+B82htDAJwN/EvJuANor1l+I3BIyZhjgO7i9r8klDW74/b7gS+Uea3SXqL3h8BvBv7Ohyjno8AHU/8ugNcXz/nU4vtdytSv4+efAs4tWdZ1wBuL28cBN5aIuRU4orj9TuDjqa//yXyMdw/6UODBiHgoIjYDVwAnjhQUg28iMFLM8oi4o7i9BlhCu9EZLiYiYm3x7ZTiKHVVVdKetFf5uyilnqmKHsfhwMUAEbE5IlYlPMSRwK8j4jclzu0Gpkvqpt3gPl4i5gDg5ohYHxF9wE+BkwaeNMTv9ETabz4UX99cJi4ilkTEfUNVaIiY64r6AdwM7FkyrnNb8hkMeH0M81r9NPChgeePEDOsIeLeA5wXEZuKc1aWLUuSgD8FLi9ZVgDP9YDnMOD1MUSMN/8Yxng30HsAj3Z8v4wRGs06SNoHOIh2j3ikc7uKj3grgesjYsSYwr/R/gNMWa0lgOsk3S6pbIbVfsCTwBeL4ZSLJM1IKPMUBvkD3K5iEY8BnwQeAZYDz0bEdSUefzFwuKSdJfXS7lkNucHDALtGxPKi/OXALiXjqnon8P2yJ0v6J0mPAm8Fzi1x/gnAYxFxd2K9ziiGUy4ZbLhnCC8BXivpFkk/lfR7CeW9FlgREQ+UPP8DwCeK5+KTwIdLxJTe/GNHNN4NtAa5b1Tn/UmaCXwL+MCA3s+gImJrRCyi3aM6VNKBJco4HlgZEbcnVu81EXEw8EbgvZIOLxHTTftj4+cj4iBgHe3hgBFJ6qH9x/HNEufOo92j3RfYHZgh6W0jxUXEEtpDBtcDP6A9jNU3bNA4knQO7fp9tWxMRJwTEXsVMWeM8Pi9wDmUaMgH+DzwImAR7TfIT5WM6wbmAYfR3jv0G0XPuIxTKfHm3eE9wFnFc3EWxae6EbyT9mv9dtpDj0Nu/rEjGu8GehnPf8fck3Ifm7NImkK7cf5qRFyVElsMG9zI9hvoDuY1wAmSHqY9bPMGSV8pUcbjxdeVtDfjPXT4CKD9HC7r6NlfSbvBLuONwB0RsaLEuUcB/xMRT0bEFtrLzL66TCERcXFEHBwRh9P+iFu2R7ZC0gKA4uvKEc6vRNJpwPHAWyMip6PwNUb+iP4i2m9ydxevjz2BOyTtNlxQRKwoOgv9wIWUe21A+/VxVTFc9wvan+jmjxBDMYz1x8DXS5YDcBrt1wW03/RHrGNELI2IYyLilbTfDH6dUN6kN94N9K3AQkn7Fr25U2hvoVW7otdwMbAkIv61ZMwLnruaL2k67UZq6UhxEfHhiNgzIvah/X/6cUQM29uUNEPSrOdu075oNeIslYh4AnhU0v7FXUcCvxoprpDSQ3oEOExSb/FcHkl7HH9EknYpvu5N+4++bJlX0/6jp/j63ZJxySQdC/w9cEJErE+IW9jx7QmM8PqIiF9GxC4RsU/x+lhG++L1EyOUs6Dj25Mo8doofAd4Q/EYL6F9EfmpEnFHAUsjYlnJcqDduTqiuP0GSrwRd7w2vPnHYMb7KiXtMcn7ab9znlMy5nLaH/O20H6Bv6tEzB/QHj65B7irOI4bIeblwJ1FzGIGuZpdotzXUWIWB+2x5LuL496yz0URuwi4rajnd4B5JWJ6gaeBOQnl/CPtBmgx8GWKmQEl4v6L9pvG3cCRZX+nwM7ADbT/0G8AdioZd1JxexOwAvhhiZgHaV8Pee618YWSZX2reD7uAb4H7JHyWmWQmTtDlPNl4JdFOVcDC0rWrwf4SlHHO4A3lKkfcCnw1yl/g7T/xm4vfs+3AK8sEXMm7b//+4HzKLKbfbQPp3qbmTXUeA9xmJnZENxAm5k1lBtoM7OGcgNtZtZQbqDNzBrKDbSZWUO5gTYza6j/D3Gg+vdLONx/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the correlation matrix.\n",
    "# vmin and vmax set the color scale.\n",
    "sns.heatmap(np.corrcoef(returns_arr), cmap=\"PiYG\", vmin=-1, vmax=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Capital Asset Pricing Model as an Example\n",
    "\n",
    "In the Capital Asset Pricing Model, we attempt to explain stocks' excess returns with the excess returns of the market. Before we go into the finer details of the assumptions, let's run a regression on the returns of Apple stock (\\$AAPL). Download the data file for this week's lecture from Blackboard and load it.\n",
    "\n",
    "Use the command \n",
    "\n",
    "```python\n",
    "df = pd.read_csv(\"./Data_Lecture_3.csv\", index_col=0, parse_dates=True)\n",
    "```\n",
    "\n",
    "to create a DataFrame. The file is more comprehensive than what we need for now. It contains the monthly returns of most of the stocks in the S\\&P 500 for the last years. In addition, it contains the returns for some sample portfolios -- like the market portfolio --, the risk-free rate and the industrial production and CPI changes. But we'll come to that later.\n",
    "\n",
    "At its core, the CAPM is a single-factor model, and we're going to treat it like that. Our equation we wish to estimate is\n",
    "\n",
    "$$R_{AAPL,t}^e = b_{AAPL} F_t + \\sigma_{AAPL} \\epsilon_{AAPL,t},$$\n",
    "\n",
    "where we take $F_t$ to be the market portfolio. When we use this approach, then we separate the risk of Apple into two categories: the market risk that we're taking by buying the stock and the risk that is specific to the stock: \n",
    "\n",
    "$$b_{AAPL}^2 \\sigma_F^2 +\\sigma_{\\epsilon,C}^2.$$\n",
    "\n",
    "Let's get to the estimation and interpretation.\n",
    "\n",
    "#### Problem 2a.\n",
    "\n",
    "In the first step, select from the DataFrame the columns with Apple's return, the risk-free rate, and the market portfolio return.\n",
    "Remember: columns are accessed and created using the ```.loc[]``` command. \n",
    "The market portfolio return is already net of the risk-free rate, but Apple's return is not. Therefore, let's create a new column \"AAPL_E\" which is the excess return.\n",
    "\n",
    "To run a regression to estimate $b_{AAPL}$, we need another module.\n",
    "If you have not already installed it, get ```statsmodels``` and import it as follows:\n",
    "\n",
    "```python \n",
    "import statsmodels.formula.api as sm\n",
    "```\n",
    "\n",
    "The syntax for the regression model is\n",
    "\n",
    "```python\n",
    "mod = sm.ols(formula=\"Y~X\", data=df).fit()\n",
    "mod.summary()\n",
    "```\n",
    "\n",
    "which first models and fits $Y_t = a + bX_t + \\epsilon_t$ with the data in ```df```, saves the results and outputs a summary table.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>AAL</th>\n",
       "      <th>AAP</th>\n",
       "      <th>AAPL</th>\n",
       "      <th>ABC</th>\n",
       "      <th>ABT</th>\n",
       "      <th>ACN</th>\n",
       "      <th>ADBE</th>\n",
       "      <th>ADI</th>\n",
       "      <th>ADM</th>\n",
       "      <th>...</th>\n",
       "      <th>XYL</th>\n",
       "      <th>YUM</th>\n",
       "      <th>ZBH</th>\n",
       "      <th>ZION</th>\n",
       "      <th>MKT</th>\n",
       "      <th>SMB</th>\n",
       "      <th>HML</th>\n",
       "      <th>RF</th>\n",
       "      <th>CPI</th>\n",
       "      <th>IndProd</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2013-01-31</th>\n",
       "      <td>8.965403</td>\n",
       "      <td>5.617027</td>\n",
       "      <td>1.604202</td>\n",
       "      <td>-15.559467</td>\n",
       "      <td>4.947367</td>\n",
       "      <td>7.786493</td>\n",
       "      <td>7.793523</td>\n",
       "      <td>0.397299</td>\n",
       "      <td>3.687699</td>\n",
       "      <td>4.077818</td>\n",
       "      <td>...</td>\n",
       "      <td>3.016765</td>\n",
       "      <td>-1.707844</td>\n",
       "      <td>11.253543</td>\n",
       "      <td>8.592044</td>\n",
       "      <td>5.57</td>\n",
       "      <td>0.39</td>\n",
       "      <td>0.95</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.281588</td>\n",
       "      <td>0.405392</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2013-02-28</th>\n",
       "      <td>-7.655023</td>\n",
       "      <td>-6.136895</td>\n",
       "      <td>3.763957</td>\n",
       "      <td>-2.577849</td>\n",
       "      <td>4.403527</td>\n",
       "      <td>-0.265997</td>\n",
       "      <td>3.378099</td>\n",
       "      <td>3.837650</td>\n",
       "      <td>4.302967</td>\n",
       "      <td>11.621178</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.140784</td>\n",
       "      <td>0.828099</td>\n",
       "      <td>0.481413</td>\n",
       "      <td>3.537714</td>\n",
       "      <td>1.29</td>\n",
       "      <td>-0.45</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.209016</td>\n",
       "      <td>-0.172034</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2013-03-31</th>\n",
       "      <td>1.459880</td>\n",
       "      <td>23.395607</td>\n",
       "      <td>8.016544</td>\n",
       "      <td>0.285049</td>\n",
       "      <td>8.621657</td>\n",
       "      <td>4.428448</td>\n",
       "      <td>2.142036</td>\n",
       "      <td>10.162677</td>\n",
       "      <td>2.769777</td>\n",
       "      <td>5.703634</td>\n",
       "      <td>...</td>\n",
       "      <td>0.217944</td>\n",
       "      <td>9.408769</td>\n",
       "      <td>0.615394</td>\n",
       "      <td>3.419136</td>\n",
       "      <td>4.03</td>\n",
       "      <td>0.78</td>\n",
       "      <td>-0.29</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.041407</td>\n",
       "      <td>0.109350</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2013-04-30</th>\n",
       "      <td>-1.270848</td>\n",
       "      <td>-0.413346</td>\n",
       "      <td>1.477238</td>\n",
       "      <td>0.027105</td>\n",
       "      <td>5.059334</td>\n",
       "      <td>4.805337</td>\n",
       "      <td>7.999121</td>\n",
       "      <td>3.533298</td>\n",
       "      <td>-5.527490</td>\n",
       "      <td>0.620661</td>\n",
       "      <td>...</td>\n",
       "      <td>0.687039</td>\n",
       "      <td>-4.955616</td>\n",
       "      <td>1.621978</td>\n",
       "      <td>-1.491662</td>\n",
       "      <td>1.55</td>\n",
       "      <td>-2.42</td>\n",
       "      <td>0.63</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.237758</td>\n",
       "      <td>0.203532</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2013-05-31</th>\n",
       "      <td>9.236622</td>\n",
       "      <td>3.887928</td>\n",
       "      <td>-2.853882</td>\n",
       "      <td>2.224061</td>\n",
       "      <td>0.311349</td>\n",
       "      <td>-0.679443</td>\n",
       "      <td>0.819326</td>\n",
       "      <td>-4.933379</td>\n",
       "      <td>5.058140</td>\n",
       "      <td>-4.613367</td>\n",
       "      <td>...</td>\n",
       "      <td>1.797858</td>\n",
       "      <td>-0.544640</td>\n",
       "      <td>2.658907</td>\n",
       "      <td>13.259207</td>\n",
       "      <td>2.80</td>\n",
       "      <td>1.67</td>\n",
       "      <td>2.60</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.195554</td>\n",
       "      <td>-0.427080</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 455 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                   A        AAL       AAP       AAPL       ABC       ABT  \\\n",
       "Date                                                                       \n",
       "2013-01-31  8.965403   5.617027  1.604202 -15.559467  4.947367  7.786493   \n",
       "2013-02-28 -7.655023  -6.136895  3.763957  -2.577849  4.403527 -0.265997   \n",
       "2013-03-31  1.459880  23.395607  8.016544   0.285049  8.621657  4.428448   \n",
       "2013-04-30 -1.270848  -0.413346  1.477238   0.027105  5.059334  4.805337   \n",
       "2013-05-31  9.236622   3.887928 -2.853882   2.224061  0.311349 -0.679443   \n",
       "\n",
       "                 ACN       ADBE       ADI        ADM  ...       XYL       YUM  \\\n",
       "Date                                                  ...                       \n",
       "2013-01-31  7.793523   0.397299  3.687699   4.077818  ...  3.016765 -1.707844   \n",
       "2013-02-28  3.378099   3.837650  4.302967  11.621178  ... -1.140784  0.828099   \n",
       "2013-03-31  2.142036  10.162677  2.769777   5.703634  ...  0.217944  9.408769   \n",
       "2013-04-30  7.999121   3.533298 -5.527490   0.620661  ...  0.687039 -4.955616   \n",
       "2013-05-31  0.819326  -4.933379  5.058140  -4.613367  ...  1.797858 -0.544640   \n",
       "\n",
       "                  ZBH       ZION   MKT   SMB   HML   RF       CPI   IndProd  \n",
       "Date                                                                         \n",
       "2013-01-31  11.253543   8.592044  5.57  0.39  0.95  0.0 -0.281588  0.405392  \n",
       "2013-02-28   0.481413   3.537714  1.29 -0.45  0.03  0.0 -0.209016 -0.172034  \n",
       "2013-03-31   0.615394   3.419136  4.03  0.78 -0.29  0.0  0.041407  0.109350  \n",
       "2013-04-30   1.621978  -1.491662  1.55 -2.42  0.63  0.0  0.237758  0.203532  \n",
       "2013-05-31   2.658907  13.259207  2.80  1.67  2.60  0.0  0.195554 -0.427080  \n",
       "\n",
       "[5 rows x 455 columns]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df=pd.read_csv('./data/Data_Lecture_3.csv',index_col=0,parse_dates=True)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alpakpinar/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:4: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  after removing the cwd from sys.path.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>AAPL</th>\n",
       "      <th>MKT</th>\n",
       "      <th>RF</th>\n",
       "      <th>AAPL_E</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2013-01-31</th>\n",
       "      <td>-15.559467</td>\n",
       "      <td>5.57</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-15.559467</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2013-02-28</th>\n",
       "      <td>-2.577849</td>\n",
       "      <td>1.29</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-2.577849</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2013-03-31</th>\n",
       "      <td>0.285049</td>\n",
       "      <td>4.03</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.285049</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2013-04-30</th>\n",
       "      <td>0.027105</td>\n",
       "      <td>1.55</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.027105</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2013-05-31</th>\n",
       "      <td>2.224061</td>\n",
       "      <td>2.80</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.224061</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 AAPL   MKT   RF     AAPL_E\n",
       "Date                                       \n",
       "2013-01-31 -15.559467  5.57  0.0 -15.559467\n",
       "2013-02-28  -2.577849  1.29  0.0  -2.577849\n",
       "2013-03-31   0.285049  4.03  0.0   0.285049\n",
       "2013-04-30   0.027105  1.55  0.0   0.027105\n",
       "2013-05-31   2.224061  2.80  0.0   2.224061"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Returns of AAPL, market portfolio return, risk-free rate\n",
    "appl_rets = df[['AAPL', 'MKT', 'RF']]\n",
    "# APPL excess returns\n",
    "appl_rets['AAPL_E'] = appl_rets['AAPL'] - appl_rets['RF']\n",
    "appl_rets.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>AAPL_E</td>      <th>  R-squared (uncentered):</th>      <td>   0.258</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared (uncentered):</th> <td>   0.249</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th>          <td>   28.20</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Fri, 07 Feb 2020</td> <th>  Prob (F-statistic):</th>          <td>9.37e-07</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>10:10:38</td>     <th>  Log-Likelihood:    </th>          <td> -271.17</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    82</td>      <th>  AIC:               </th>          <td>   544.3</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    81</td>      <th>  BIC:               </th>          <td>   546.7</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>              <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>              <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "   <td></td>      <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>MKT</th> <td>    1.1060</td> <td>    0.208</td> <td>    5.310</td> <td> 0.000</td> <td>    0.692</td> <td>    1.520</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>16.283</td> <th>  Durbin-Watson:     </th> <td>   2.075</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td> <th>  Jarque-Bera (JB):  </th> <td>  21.381</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.899</td> <th>  Prob(JB):          </th> <td>2.28e-05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 4.739</td> <th>  Cond. No.          </th> <td>    1.00</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                                 OLS Regression Results                                \n",
       "=======================================================================================\n",
       "Dep. Variable:                 AAPL_E   R-squared (uncentered):                   0.258\n",
       "Model:                            OLS   Adj. R-squared (uncentered):              0.249\n",
       "Method:                 Least Squares   F-statistic:                              28.20\n",
       "Date:                Fri, 07 Feb 2020   Prob (F-statistic):                    9.37e-07\n",
       "Time:                        10:10:38   Log-Likelihood:                         -271.17\n",
       "No. Observations:                  82   AIC:                                      544.3\n",
       "Df Residuals:                      81   BIC:                                      546.7\n",
       "Df Model:                           1                                                  \n",
       "Covariance Type:            nonrobust                                                  \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "MKT            1.1060      0.208      5.310      0.000       0.692       1.520\n",
       "==============================================================================\n",
       "Omnibus:                       16.283   Durbin-Watson:                   2.075\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               21.381\n",
       "Skew:                          -0.899   Prob(JB):                     2.28e-05\n",
       "Kurtosis:                       4.739   Cond. No.                         1.00\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Regression\n",
    "import statsmodels.api as sm\n",
    "\n",
    "mod = sm.OLS(appl_rets['AAPL_E'], appl_rets['MKT']).fit()\n",
    "mod.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to know find out a couple of things about the behavior of Apple stock. \n",
    "<ul>\n",
    " <li>How does the excess return of Apple compare to the market? Does it fluctuate more or less?\n",
    " <li>Historically, the annual excess return of the market, or market risk premium, has been $\\lambda_F=\\lambda_{MKT}=6\\%$. What annual return do you expect from Apple stock?\n",
    " <li>We know that the risk is separated into two parts. Using the \"MKT\" column in the DataFrame, find an estimate for $\\sigma_F^2=\\sigma_{MKT}^2$.\n",
    " <li>The second part is the error from the residuals, $\\sigma_{\\epsilon,AAPL}^2$, which can be obtained from the estimated model through ```mod.mse_resid```.\n",
    "</ul>\n",
    "\n",
    "To learn more about the errors, make a scatter plot of the excess market returns on the x-axis and Apple's excess returns on the y-axis. The general syntax is ```df.plot.scatter(x=X, y=Y)```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'AAPL_E')"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZQddZ338fc3TYuNIo0kStIQCB6MsiiRPtGZ6IyAGEQOtCAI43HwUSfqI+dRn8dIIi4geBLNOO6KURlxdARkaQJBgywOo6NoZ4EAIRplSydiswQwtBCS7/NHVadv3666a2333s/rnD7crqp775fqSn3rt5u7IyIiUm5K3gGIiEgxKUGIiEgkJQgREYmkBCEiIpGUIEREJJIShIiIRMo1QZjZgWZ2q5ltMLO7zezD4fYXm9nPzewP4X/3zTNOEZFOZHmOgzCz6cB0d19jZnsDq4EB4N3AY+6+1MwWAfu6+7m5BSoi0oFyLUG4+1Z3XxO+fgrYAPQBpwCXhoddSpA0REQkQ7mWIEqZ2cHAbcARwIPu3luy73F3r1jNNHXqVD/44IPTDFFEpO2sXr36EXefFrVvj6yDiWJmLwSuAj7i7k+aWa3vWwAsAJg5cyZDQ0PpBSki0obM7IG4fbn3YjKzboLk8CN3vzrc/HDYPjHWTvGXqPe6+3J373f3/mnTIhOgiIg0KO9eTAZ8D9jg7v9WsmsFcHb4+mzg2qxjExHpdHlXMc0D3gWsN7N14bZPAEuBK8zsvcCDwOk5xSci0rFyTRDu/ksgrsHhuCxjERGRiXJvgxARkWLKu4pJRDrI4Nphlq3ayJZto8zo7WHh/NkMzOnLOyyJoQQhIpkYXDvM4qvXM7pjJwDD20ZZfPV6ACWJglIVk4hkYtmqjbuTw5jRHTtZtmpjThFJNUoQIpKJLdtG69ou+VOCEJFMzOjtqWu75E8JQkQysXD+bHq6uyZs6+nuYuH82TlFJNWokVpEMjHWEK1eTK1DCUJEMjMwp08JoYWoiklERCIpQYiISCQlCBERiaQEISIikZQgREQkkhKEiIhEUoIQEZFIShAiIhIp9wRhZpeY2V/M7K6Sbeeb2bCZrQt/TswzRhGRTpR7ggC+D5wQsf1L7n5U+HNDxjGJiHS83BOEu98GPJZ3HCIiMlHuCaKCc8zszrAKat+oA8xsgZkNmdnQyMhI1vGJiLS1oiaIbwEvA44CtgJfjDrI3Ze7e7+790+bNi3L+ERE2l4hE4S7P+zuO919F/AdYG7eMYmIdJpCJggzm17y69uAu+KOFRGRdOS+HoSZ/Rh4IzDVzDYDnwHeaGZHAQ7cD7w/twBFRDpU7gnC3c+K2Py9zAMREZEJClnFJCIi+cu9BCEikqfBtcNaJzuGEoSIdKzBtcMsvno9ozt2AjC8bZTFV68HUJJAVUwi0sGWrdq4OzmMGd2xk2WrNuYUUbGoBCEidWu2WqYo1Tpbto3Wtb3TKEGISF2arZYpUrXOjN4ehiOSwYzenkzjKCpVMYlIXZqtlilStc7C+bPp6e6asK2nu4uF82c3/dmDa4eZt/QWZi1aybyltzC4drjpz8yaShAiUpdmq2WKVK0zVmJJurprcO0wC6+8gx07HQhKSQuvvGPCd7YCJQgRqUuz1TJFq9YZmNOX+E37guvu3p0cxuzY6Vxw3d0tlSBUxSQidWm2WibNap2iePzpHXVtLyqVIESkLs1Wy6RVrSPJU4IQkbo1Wy2TRrVOkfT2dLNtdHJpobenO4doGqcqJhGRhJ1/8uF0T7EJ27qnGOeffHhOETVGJQgRkYS1SzWaEoSISAraoRpNCUJEdivKFBhSDEoQIgIUawoMKYbcE4SZXQKcBPzF3Y8It70YuBw4mGDJ0TPc/fG8YhTJSp5P8JWmwBiLQSWMzlKEXkzfB04o27YIuNndDwVuDn8XaWtjT/DD20Zxxp/gs5rDp9oUGHnHJ9nLPUG4+23AY2WbTwEuDV9fCgxkGpRIDvKexC5uqoux7XnHl5d2mHSvUblXMcV4qbtvBXD3rWb2kqiDzGwBsABg5syZGYYn0pyoqpq8J7FbOH/2hDYImDgFRt7x5aHT22VyL0E0w92Xu3u/u/dPmzYt73BEahJXVdO7V/Qo26wmsRuY08eSU4+kr7cHA/p6e1hy6pG7b4TVShhpyPvpvVNLTWOKWoJ42Mymh6WH6cBf8g5IJClxN50995hCT3dX7BN8LZptRC7tuz/2WR+9fB0zens45hXTuGr1cFPx1aMIT++dWGoqVdQSxArg7PD12cC1OcYikqi4m8sTozsqPsFXk2QjctRnXbV6mNOO7ms4vnoV4ek9j1JTTZ58Et73PjALfn75y1S+JvcShJn9GHgjMNXMNgOfAZYCV5jZe4EHgdPzi1AkWZXWQ2hm9G0t3VSb/axb7x3hV4uObSi+WpSWgDzmmCyf3qu1y2Tqr3+Fj38cvvWtyfuOPDKVr8y9BOHuZ7n7dHfvdvcD3P177v6oux/n7oeG/y3v5STSstJaDyHJ6pA8qlbKSy1xsnx6r9Yuk7qnn4aPfCQoJey998Tk8OEPw/bt4A777JPK1+deghDpNGlN5JbkSm15rPoWVWopl8fTe+ZzKm3bBu96F1x//eR9H/wgfOEL8MIXZhKKEoS0jVYa5ZvGTSfJ6pA8qlYqlU4MCv83bcpTT8GLXhS9733vgy9+MX5/ipQgpC0UocdL3pIsmeQxXXVcqaWvtyfVdo/cPP10UG20a1fk7ld9+DJ27L1PUKWVQ3IAMPdKtX2to7+/34eGhvIOQ3Iyb+ktnXVzaUPlSR6CUkumdf5pe+YZmDYtKDFEOOZfvs19L574/5r2NWxmq929P2qfShDSFjq9v3pS8qymiyu1QPAA0ApVh5F27ICDDoKtW6P3r18PRxzBrEUrIxvn87yGlSCkLeTRqNpuilBNV942U4SYGvLss7DnnvH716yBOXMmbCriNZx7N1eRJKTVdbQVNTo9RREGppUrYkyxnntufOBaVHL49a+DLqnuk5IDFPMaVglC2kK7rAHcrLgn7qEHHuPWe0cmVd2Unq+op1eYWMWRdRVU4asOd+2Crq74/RdeCJ/8ZE0fVcRrWI3UIm0krrHeYEL9dneXgcOOXR57zJixRtJ6GpGTSiSF7HzgDlMqVL4sXBiMVWgRlRqpVcUkhZH3zJ2tbnDtcGwpoPzGv2OnT0gOY8dY2XGlVRy1VvckOSdUoapdxqqPopLDggXj1UctlByqURWTFELLNkampN4n8LHz1ywneDqP+t5aq3uSnBMq92oXK0+ZJc44Ay6/PJs4cqIEIYWQ5E2l1TWSLCtNUxFXdRSlUtVNrb1skm43yHyqi0pJ4fjj4cYbs4slZ6pikkKo56bS7lVRjfTcqXTzfefrZk6qpunuMrqnTLwRVqu6qbW6p7BTZFcyVn0UlRzMxquPOig5gBKEFEStN5Uk67eLqpEn8Ljz19fbw0UDR06akXTZ21/NstNfXdcspbXObFqodoNKKiUFGE8KMVNhdAJVMUkh1Do5XCdURTUyYKra+YurpmmkTaDae7JsN6i7t1Sl6iMIEoLspgQhhVDrTaXw/eIT0MhMqrk35kbEk/Z319xWo6TQMCUIKYxabipFnI4gaY3e7DNvzE1BPSWCiqXJ1xxQ+YvaJCmkPXCx0AnCzO4HngJ2As/FDeaQzlGoJSBT1A43+3rV23urvNR4/+dPqvwFbZIUxmTRNbzQCSJ0jLs/kncQUgxFq0pJWistepS0etuXZvT28KvFx1X+0DZLCqWyaI9rhQQhMkG7Pl3nMVhwLCENbxuly4yd7vTllJhqbl8K2xR+FfM5g6sfql7F1AayaI8rejdXB240s9VmtqB8p5ktMLMhMxsaGRnJITyR5KQxc2mlMSOlXYYBdoZP23l1Ha7Y1blKl9RDPr6CeUtuZnDN5o5IDpDNeJOilyDmufsWM3sJ8HMzu9fdbxvb6e7LgeUQTNaXV5AiSUj6ibBaiaTS6Os0ug5Xqz4rb1+q2qbw7LPQ3Q3AnxKLsnVk0R5X6ATh7lvC//7FzK4B5gK3VX6XSGtKuodWtTrqaoknyaqKWqrPBub0VX/6374d9torsbhaWRbtcYVNEGb2AmCKuz8Vvn4z8NmcwxJJTdJPhNVKJJXWgBjbn5SmuqSOjMDUqYnF0k7Sbo8rbIIAXgpcY0Gd4x7Af7r7z/INSSQ9ST8RViuRRCWkMUlXVdTdJXXjRnj5yxP7fmlMYROEu/8JeHXecYhkKcknwlqm3wAy6cVUU5fU22+HuXMT+05pXmEThIg0p5YSSepdhqdMAffYLqn/87Uf8PfnvCu975emKEGItLFaEkDig/P6+mDLltjdn3nT+7npuDM6ahBgq1KCkI7VyaOWxyQ2OO+YY+AXv4jfv2gRLFkCwAXhjxRf0QfKiaSiE9aVqEUjg/PGBt+tfMUbxgevRSSHrf/4ZuYtuZlZ517PvH2O77hz2w5UgpCO1AnrStSi3sF5D558BgPX/YSBuA985SvhnnvGSybh53T6GuOtSglCOlKrrSuRVnVYTYPzLroIPvUpAGbGfM7LFq1kl3sQWxirEnDrU4KQjtRK60qkOYlfXFfYL/sGsMrdUg8+9/rxX8rmcYqbwqOoCViiKUFIR2qldSUuuO7u1J7GS7vCzrzjdn582ScqHj9vyc0VR1+PxWYWPdN2EROwxFOCkI7UKutKDK4d5vGnd0TuS+Rp/M47GXjNq+PbFGDCnX5hWWmmhrfs1t1lhUzAEk8JQjpWK6wrUak3UcNP4w8+CAcdVPmYmIV2yhPrlHD0dS1e8Lw9Cn++ZaKK3VzN7IqS158v23djWkGJSKBSKaGup/Ft28a7pMYlB/fxnwoG5vTxq0XHct/St/LFM15NT3dXTSE8MRpdEpLiqjYO4tCS18eX7ZuWcCwiUiaulNDb0139afzZZ8eTwr77Rh9TY1KIMzCnjyWnHklfbw8G9PX2sO9e3ZHHqv2h9VRLEJWuGi3QI5KyhfNnT3pC7+nu4vyTD49+g/t4Uthzz/hjmkgK1bz1VdMjY1b7Q+up1gaxl5nNIUgkPeFrC3/0OCCSspob02OW4twtpWQQ1QX3qtXDnHZ0H7feO1LoDgBSnXmFC8fMbq30Znc/JvGIGtTf3+9DQ0N5hyGSmvLBclWnz04pKZSat/SWyG6vfb09/GrRsal/vzTPzFa7e3/UvooliFoTgJkd7+4/byQ4kby00mR9Y0/qGy56S+UDM0gKpdIakd5Kf5t2llQ3188DiScIMzsB+ArQBXzX3Zcm/R3SmdIcnZw4MwYgfqxCxkmhVBoj0lvqb9PmkkoQVSpAG/hAsy7gGwS9pzYDvzOzFe5+T9LfJZ2n8HMFVWlTGJvmwoD7SrZXe/JO+sk8jRHphf/bdJCkEkQajzBzgU3h0qOY2WXAKYAShDStKJP1ld6w76uyTvOEuY9CpU/q1Z6803gyT2NEelH+NlLskdR9wEMlv28GXptTLNJmijBZ3+DaYQZec0BN01wMrh2mp8qTerUn77SezJMekV6Ev40EGl4wyMxKb9b3Nx/K5K+I2DahpGJmC8xsyMyGRkZGUghB2lXc+IJM+uqH4xQGXnNA5O55S26eNE5hYE4fpx3dR1dY9dRlxmlHT7wxV3vybpUn81z/NjJBMyWInxBOD+/upyYTzgSbgQNLfj8AmLDQrbsvB5ZD0M01hRikTWU+WV+NbQoAFnHDHlw7zFWrh3fPe7TTnatWD9N/0It3x1ztybtVnsxbZSLFTtBMgki8YbrM74BDzWwWMAycCfxTyt8pHST1yfqqJIW4qbOnmAXVT2Fsg2uH+X9X3DFpUrzy6qFqDcatNMV5K0yk2AmaSRCpPrG7+3Nmdg6wiqCb6yXufnea3ynStGojmnft2n1M3NTZO913Nx4DLL56feyMqaXVQ9WevPVkPpnGW1RWbST1dUQnAgOOdfcXpBVYvTSSWnJTLSns2AF7RD+LxZUOIBiNDFRcoEcjlhtX3qsLghLVklOP7Kgk0fBIauBfG9wn0t6qJYXt22Gvvap+zMCcPj56+brIfdVWbiuvHtLTcH003qK6alNt/Fc4Qd/LgLvdfUM2YYk0p9abZV031WpJ4ZFHYL/96o41rvG4ki6zCU+6Gn1cv1bp1ZWnagsGfRq4HDgNWGlm/5JJVCJNGLtZDm8bxRm/WQ6uHa7/uLGps+OSwwMPjHdJbSA5QHS3zkp6urv44hmvnnDjr/Q0LNHiem8VrVdXnqpVMb0DOMrdnzaz/YCfAd9JPyyReNWe+mutOog7Lm58wm733AOvfGXz/yOh8sbjSr0/+mJKOa38NJxk1Vg9n9VKvbryUi1B/M3dnwZw90fNrOGBdSJJqKUqpdabZenv/33xeznwiYfjv/jXv4bXva6Z0Csq7dbZyBTarTLGoVySVWP1fpZ6dVVXLUG8zMxWhK+t7Hfc/eTUIhOJUEvpoNab5bLbvsvbfz0Y/2WrVsGb39x80HVq5Mk26j3dXcb2Z55j1qKVhb35JdlQ3MhnabxFZdUSxCllv6vnkuSqltJBxRvsRRfBpz4FwNsjPudjAx/n9Z/+P7neNBp5si1/T+9e3fz1b8+xbXQHUNxG6ySrxlq5mq2oqvZiitpuZgcSjGyO3C+SllpKB+U3yw/8/mbOveZLcFH0Z37u7Qv57sv+sVBP2Y082ZZXUz3+9I4J+4vYhTPJqrFWrWYrsppHUpvZVOB04CyCmVavSSsokTi1Vr8MbFnHwOIK02d//evwoQ8BcF74U4+ijzlolafpJBuK1eicvIoJwsz2Bt5GMAfSywmSwiHuXqWbh9Si6DeZIqpY/fK738HcufFv/vSn4YILmo6hFcYcxD1NTzErVJtEkg3FanROXrWpNkaB3wKfBH7p7m5mf3L3Q7IKsFatNtWGhvkn5N57K3c5ff/74eKLE/3KRnoZZS3q+iqn602g8lQb1bqtfgJ4PvAtYLGZvSzp4DqVBjY1YcuW8cFrUcnhox8dH7zWRHIYXDvMvKW3MGvRSuYtvWX3ALpWqL4ZmNPHklOPpK+3B4Pd60iU0vUm1VRrpP4S8CUzO4Sg7WEQmGFm5wLXuPvvM4ixLbXCTaZQHn0Upk6N3//Od8IPf5jY11WqRmqVxtDSRutZi1ZGHqPrTSqpaeCbu//J3T/n7kcSrBU9H/hpqpG1OQ3zr8H27eMlhajkcNxx4yWFBJMDVC7hteKKZ7repBE1j4w2s6PM7PPAivB9X04tqg7QijeZTDz77HhSeOELJ+8//vjxpHDTTamFUamEV15909fb03Bdflw1VtJ0vUkjqvViejnBeIezgEcJJu4zd39j+qG1N/W4KLFzZ+x6CQCcdBJcd1128VC9GimJEbhZ9obS9SaNqNaLaRfw38B73X1TuE29mKR57jClQgF27ly4/fbs4imTRS+zVugNJe2vmV5MpwF/Bm41s++Y2XGkvxY1Zna+mQ2b2brw58S0v1My4D5efRSVHEqrj3JMDjC5F1Az1Uhx1FFBiq5aL6ZrgGvM7AXAAPBR4KVm9i2CXkw3phjbl9xdcz+1g0oL7cyZA2vWZBdLHdKeyK1VekNJ56q1F9N2d/+Ru58EHACsAxalGpm0tkoL7Rx44HhJoaDJIQtqOJaiq3t9B3d/zN2/7e5pV5KeY2Z3mtklZrZvyt8lSaiUFPbcczwpPPhg9rEVUBbVWCLNqNhIneoXm90E7B+x6zzgN8AjgAMXAtPd/T0Rn7EAWAAwc+bMox944IH0ApZo1dZpzun6EpHaVGqkzi1B1MrMDgaud/cjKh2nXkwZUlLInSZ6lKRUShA1T/edJTOb7u5bw1/fBtyVZzxCSyeFdruZtsJsstIeCpkggC+Y2VEEVUz3A+/PN5wO1cJJYUw73kyTXKazXu2WbKWyQiYId39X3jF0rGnT4JFH4ve3QFIolefNNC15jZ9ox2QrldXdi0na0GGHjfc+ikoOY72PWiw5QHsORstr4j1NUd95lCA6SOnEcCvmvnU8KWzYMPngFk4KpdpxFtO8xk+0Y7KVygpZxSTJG1w7zB3nLeFXP/1G/EG7dlVvd2gx7bhOcV4T72nkd+dRgmh3P/gBnH02AwRzpZR7/ed+zi8/8aaso8pMu85imvY0IFHaMdlKZUoQ7ejKK+H002N3H7LwWnZNCaoo7MlnsooqN3ncTGvVSr2C2jXZSjwliHZxww3w1rfG7v6Hi1bx4FM7Jm1X9UB+WrFXUJGTrSRPjdSt7NZbxxuao5LD3/62u6H5/554eGTD5jGvmJbJimYymXoFSdGpBNFqNm6EV7wifv/27bDXXpM2R1UPHPOKaVy1erilnmDbiXoFSdEpQbSC++6DQyos4vfEE/CiF1X9mPLqgXlLb2m7QWStRL2CpOhUxVRUmzePVx9FJYdHHx0fp1BDcoiiJ9h8aT0IKTqVIIrkz3+G6dPj92/bBvvsk9jX6Qk2X+oVJEWnBJG3v/4V9t47fv8jj8B++6Xy1erXnj/1CpIiUxVTHv72N/j4x4Pqo6jk8Oc/j1cfpZQcQCuaiUhlKkFk5Zln4MIL4XOfi97/0ENwwAHZxoSeYEUknkoQadqxAy64ICgpPP/5E5PD2WfD44+PlxRySA4iIpWoBJG0556DZcvgE5+YvO+ss+BrX0u12khEJClKEEnYuRO+/GX42Mcm7zvtNPjmN+ElL8k+LhGRJuRWxWRmp5vZ3Wa2y8z6y/YtNrNNZrbRzObnFWNFu3bBV78aVB/tscfE5HDyybB1a1B1dOWVSg4i0pLyLEHcBZwKfLt0o5kdBpwJHA7MAG4ys5e7+87JH5Exd/j2t+GDH5y874QT4DvfUVuCiLSN3BKEu28AsMkL1JwCXObuzwD3mdkmYC7w62wjDLnDv/87vPe9k/cdeyxccgkcdFD2cYmIpKyIvZj6gIdKft8cbpvEzBaY2ZCZDY2MjCQXgTv88IdB9dGUKROTwxveAJs2BcfcfLOSg4i0rVRLEGZ2E7B/xK7z3P3auLdFbItcGNndlwPLAfr7+5tfPPmKK+Ad75i8/bWvhUsvhdkaYSwinSPVBOHujaxluRk4sOT3A4AtyUQUwT0oBTz00MTtc+YEpYjDDkvtq0VEiqyIVUwrgDPNbE8zmwUcCvw2tW977LHx5HD44bBuXZA01qxRchCRjpZbI7WZvQ34GjANWGlm69x9vrvfbWZXAPcAzwEfSrUH0377BQlBREQmyLMX0zXANTH7PgfETFoknWZw7bCmxBbJgUZSS6ENrh2eMCV5EsuiKuGI1KaIbRAiuy1btTF2WdRGjCWc4W2jOOMJZ3DtcALRirQXJQgptKSXRU064Yi0MyUIKbS45U8bXRZV63CL1E4JQgpt4fzZ9HR3TdjWzLKoSScckXbW8Y3UarAstrG/RVJ/oyTW4dY1I52ioxNEGj1kJHlJLovabMLRNSOdpKMTRKUGS/1jb1/NJBxdM9JJOroNQg2WUi9dM9JJOroEMaO3h+GIf9hFbrBU/Xe+WvGaEWlUR5cgku4hkzYN8spfq10zIs3o6AQxMKePJaceSV9vDwb09faw5NQjC/tErkFe+Wu1a0akGR1dxQTJ9pBJm+q/i6GVrhmRZnR0CaLVaJCXiGRJCaKFqP5bRLLU8VVMrSTpUcUiIpUoQbQY1X+LSFZyq2Iys9PN7G4z22Vm/SXbDzazUTNbF/5cnFeMIiKdLM8SxF3AqcC3I/b90d2PyjgekZakwZOSljzXpN4AYGZ5hSDS8jR5oKSpqL2YZpnZWjP7LzN7Q9xBZrbAzIbMbGhkZCTL+EQKQYMnmzO4dph5S29h1qKVzFt6i2YlKJNqCcLMbgL2j9h1nrtfG/O2rcBMd3/UzI4GBs3scHd/svxAd18OLAfo7+/3pOIWaRUaPNk4lb6qSzVBuPubGnjPM8Az4evVZvZH4OXAUMLhibQ8TR7YOE3dXl3hqpjMbJqZdYWvDwEOBf6Ub1QixaTBk41T6au6PLu5vs3MNgN/B6w0s1Xhrn8A7jSzO4ArgQ+4+2N5xSlSZJo8sHGauqY6c2+Pqvv+/n4fGlItlIjUprwNAoLSV6clWDNb7e79Ufs0klpEOpKmrqlOCUJEOpamrqmscI3UIiJSDEoQIiISSQlCREQiKUGIiEgkJQgREYmkBCEiIpGUIEREJJIShIiIRFKCEBGRSEoQIiISSQlCREQiKUGIiEgkJQgREYmkBCEiIpGUIEREJFKeS44uM7N7zexOM7vGzHpL9i02s01mttHM5ucVo4hIJ8uzBPFz4Ah3fxXwe2AxgJkdBpwJHA6cAHzTzLpiP0VE6ja4dph5S29h1qKVzFt6C4Nrh/MOSQootwTh7je6+3Phr78BDghfnwJc5u7PuPt9wCZgbh4xirSjsbWYh7eN4sDwtlEWX71eSUImKUobxHuAn4av+4CHSvZtDrdNYmYLzGzIzIZGRkZSDlGkPSxbtZHRHTsnbBvdsZNlqzbmFJEUVaprUpvZTcD+EbvOc/drw2POA54DfjT2tojjPerz3X05sBygv78/8hgRmWjLttG6tkvnSjVBuPubKu03s7OBk4Dj3H3sBr8ZOLDksAOALelEKNJ5ZvT2MByRDGb09uQQjRRZnr2YTgDOBU5296dLdq0AzjSzPc1sFnAo8Ns8YhRpRwvnz6ane2K/j57uLhbOn51TRFJUqZYgqvg6sCfwczMD+I27f8Dd7zazK4B7CKqePuTuOyt8jojUYWBO0KS3bNVGtmwbZUZvDwvnz969XWSMjdfstLb+/n4fGhrKOwwRkZZiZqvdvT9qX1F6MYmISMEoQYiISCQlCBERiaQEISIikZQgREQkUtv0YjKzEeCBCodMBR7JKJxmtEqc0DqxKs7ktUqsirO6g9x9WtSOtkkQ1ZjZUFxXriJplTihdWJVnMlrlVgVZ3NUxSQiIpGUIEREJFInJYjleQdQo1aJE1onVsWZvFaJVXE2oWPaIEREpD6dVIIQEZE6tFWCMLPTzexuM9tlZv1l+xab2SYz22hm82PeP8vMbjezP5jZ5Wb2vAxivtzM1oU/95vZupjj7jez9eFxucxKaGbnm9lwSbwnxhx3QkH7O6cAAAWwSURBVHieN5nZohziXGZm95rZnWZ2jZn1xhyXyzmtdn7Cqe4vD/ffbmYHZxVbSQwHmtmtZrYh/Df14Yhj3mhmT5RcD5/OOs6SWCr+LS3w1fCc3mlmr8khxtkl52qdmT1pZh8pO6Yw5xQAd2+bH+CVwGzgF0B/yfbDgDsIphefBfwR6Ip4/xXAmeHri4EPZhz/F4FPx+y7H5ia8/k9H/hYlWO6wvN7CPC88LwflnGcbwb2CF9/Hvh8Uc5pLecH+N/AxeHrM4HLc/hbTwdeE77eG/h9RJxvBK7POrZG/pbAiQTLGhvwOuD2nOPtAv5MMAahkOfU3durBOHuG9w9amHdU4DL3P0Zd78P2ATMLT3AgkUpjgWuDDddCgykGW/E958B/Dir70zJXGCTu//J3Z8FLiM4/5lx9xvd/bnw198QrEpYFLWcn1MIrj8IrsfjwusjM+6+1d3XhK+fAjYQszZ8izgF+IEHfgP0mtn0HOM5Dviju1ca3Ju7tkoQFfQBD5X8vpnJF/t+wLaSG0vUMWl6A/Cwu/8hZr8DN5rZajNbkGFc5c4Ji+iXmNm+EftrOddZeg/Bk2OUPM5pLedn9zHh9fgEwfWZi7CKaw5we8TuvzOzO8zsp2Z2eKaBTVTtb1m06/JM4h8Gi3JOc11RriFmdhOwf8Su89z92ri3RWwr775VyzENqTHms6hcepjn7lvM7CUEq/Dd6+63JRFfrbEC3wIuJDgvFxJUib2n/CMi3pt4V7lazqmZnUewKuGPYj4mk3NaJtdrsV5m9kLgKuAj7v5k2e41BFUkfw3bowYJlgjOQ7W/ZZHO6fOAk4HFEbuLdE5bL0G4+5saeNtm4MCS3w8AtpQd8whBsXOP8Kkt6piGVIvZzPYATgWOrvAZW8L//sXMriGoqkj8Zlbr+TWz7wDXR+yq5Vw3rYZzejZwEnCch5W7EZ+RyTktU8v5GTtmc3ht7AM8lnJck5hZN0Fy+JG7X12+vzRhuPsNZvZNM5vq7pnPKVTD3zKT67JGbwHWuPvD5TuKdE6hc6qYVgBnhr1DZhFk5N+WHhDeRG4F3h5uOhuIK5Ek7U3Ave6+OWqnmb3AzPYee03QCHtXRrGVxlFaZ/u2mBh+BxxqQY+w5xEUpVdkEd8YMzsBOBc42d2fjjkmr3Nay/lZQXD9QXA93hKX5NIStnl8D9jg7v8Wc8z+Y20jZjaX4H7yaHZR7o6jlr/lCuCfw95MrwOecPetGYc6Jra2oCjndLe8W8mT/CG4aW0GngEeBlaV7DuPoPfIRuAtJdtvAGaErw8hSBybgJ8Ae2YU9/eBD5RtmwHcUBLXHeHP3QTVKHmc3/8A1gN3EvyDm14ea/j7iQS9Xv6YR6zh3+8hYF34c3F5nHme06jzA3yWIKEBPD+8/jaF1+MhOZzD1xNUwdxZch5PBD4wdq0C54Tn7g6CzgB/n9N1Gfm3LIvVgG+E53w9Jb0cM451L4Ib/j4l2wp3Tsd+NJJaREQidUoVk4iI1EkJQkREIilBiIhIJCUIERGJpAQhIiKRlCBEmmBmbmb/UfL7HmY2YmbXh7+/28y+Hr6eYmaXhtOU3B7O1vlgePzY7J0H5/N/IjJZy42kFimY7cARZtbj7qPA8cBw+UHh4KeLgW7gf7n7rnD7uwn65J+TXcgitVEJQqR5PwXeGr6OGyX7FYIJ9/55LDmIFJ0ShEjzLiOYyuX5wKuYPOvpPxHMs3Wmj88WLFJ4ShAiTXL3O4GDCUoPN0QcsgY4iLI1SESKTglCJBkrgH8lunrpXoLFoC7Pe35/kXooQYgk4xLgs+6+Pmqnu/8PwaRsK81sZqaRiTRIvZhEEuDBVO1fqXLM9WY2DfiZmb3B3fObxlmkBprNVUREIqmKSUREIilBiIhIJCUIERGJpAQhIiKRlCBERCSSEoSIiERSghARkUhKECIiEun/AxkKBIXSOJaHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,1)\n",
    "ax.scatter(appl_rets['MKT'], appl_rets['AAPL_E'])\n",
    "def bestfit(x, slope):\n",
    "    return slope*x\n",
    "bestfit_line = bestfit(appl_rets['MKT'], 1.106)\n",
    "ax.plot(appl_rets['MKT'], bestfit_line, 'r')\n",
    "ax.set_xlabel('MKT')\n",
    "ax.set_ylabel('AAPL_E')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "source": [
    "Most notably, the CAPM predicts that there is no intercept. Look at your result and check the intercept. Thinking about statistical significance, discuss!\n",
    "\n",
    "Additionally, for the implications of the CAPM about the factor loading $b$ to make sense for individual stocks, we need to hold a well diversified portfolio. Speaking of which, you should have noticed above that a large part of the risk in holding Apple stock is not market risk, but firm-specific risk. This is not good because an investor who is diversified would have less risk when they add Apple to their portfolio, and therefore, they would be willing to pay more for the stock.\n",
    "\n",
    "#### Problem 2b.\n",
    "\n",
    "In that spirit, let's look at the CAPM for a bigger portfolio! Choose the first 80 stocks of the data set, using ```.iloc[:, 0:80]``` (all rows, first 80 columns), and build an equal weight portfolio. That means, each weight is $w_i=1/80$.\n",
    "\n",
    "Find the excess return of that portfolio, and run it against the excess market return as in 2a."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alpakpinar/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:2: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  \n",
      "/Users/alpakpinar/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:3: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n",
      "/Users/alpakpinar/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:4: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  after removing the cwd from sys.path.\n",
      "/Users/alpakpinar/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:5: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  \"\"\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>Portfolio_E</td>   <th>  R-squared:         </th> <td>   0.946</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.945</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   1402.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Fri, 07 Feb 2020</td> <th>  Prob (F-statistic):</th> <td>1.76e-52</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>10:26:24</td>     <th>  Log-Likelihood:    </th> <td> -104.69</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    82</td>      <th>  AIC:               </th> <td>   213.4</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    80</td>      <th>  BIC:               </th> <td>   218.2</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -0.0282</td> <td>    0.102</td> <td>   -0.277</td> <td> 0.782</td> <td>   -0.231</td> <td>    0.175</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>MKT</th>       <td>    1.0833</td> <td>    0.029</td> <td>   37.449</td> <td> 0.000</td> <td>    1.026</td> <td>    1.141</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 1.087</td> <th>  Durbin-Watson:     </th> <td>   2.359</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.581</td> <th>  Jarque-Bera (JB):  </th> <td>   0.546</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.106</td> <th>  Prob(JB):          </th> <td>   0.761</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 3.338</td> <th>  Cond. No.          </th> <td>    3.74</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:            Portfolio_E   R-squared:                       0.946\n",
       "Model:                            OLS   Adj. R-squared:                  0.945\n",
       "Method:                 Least Squares   F-statistic:                     1402.\n",
       "Date:                Fri, 07 Feb 2020   Prob (F-statistic):           1.76e-52\n",
       "Time:                        10:26:24   Log-Likelihood:                -104.69\n",
       "No. Observations:                  82   AIC:                             213.4\n",
       "Df Residuals:                      80   BIC:                             218.2\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -0.0282      0.102     -0.277      0.782      -0.231       0.175\n",
       "MKT            1.0833      0.029     37.449      0.000       1.026       1.141\n",
       "==============================================================================\n",
       "Omnibus:                        1.087   Durbin-Watson:                   2.359\n",
       "Prob(Omnibus):                  0.581   Jarque-Bera (JB):                0.546\n",
       "Skew:                          -0.106   Prob(JB):                        0.761\n",
       "Kurtosis:                       3.338   Cond. No.                         3.74\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_lim = df.iloc[:, 0:80]\n",
    "df_lim['RF'] = df['RF']\n",
    "df_lim['MKT'] = df['MKT']\n",
    "df_lim['Portfolio'] = df_lim.apply(np.mean, axis=1)\n",
    "df_lim['Portfolio_E'] = df_lim['Portfolio'] - df_lim['RF']\n",
    "\n",
    "import statsmodels.formula.api as sm\n",
    "mod = sm.ols(\"Portfolio_E ~ MKT\", data=df_lim).fit()\n",
    "mod.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Macroeconomic Factor Models\n",
    "\n",
    "The Capital Asset Pricing Model is arguably the most important factor model, but it suffers from imperfections. An alternative are multifactor models based on macroeconomic data. If the stock price changes, then maybe this can be explained in an unexpected surprise to macroeconomic conditions to which the firm is sensitive. One simple such model is \n",
    "\n",
    "$$R_{i,t}^e = b_{i,CPI} F_{CPI,t} + b_{i,IP} F_{IP,t},$$\n",
    "\n",
    "where we use the comsumer price index and the industrial production surprises. To get there, we first have to find out those surprises.\n",
    "\n",
    "Unlike stock returns, these time series have a strong degree of autocorrelation, which means that $CPI_{t-1}$ is very useful in prediction $CPI_t$. Again, we use regression to estimate this relationship. This particular type of model is called AR(1), an autoregression with lag 1. (In reality, we'd have to build a much more sophisticated prediction with more than one lag, for example. But this is good enough for our demonstration.)\n",
    "\n",
    "#### Problem 3a.\n",
    "Run the regressions $CPI_t = \\beta_{0,CPI} + \\beta_{1,CPI} CPI_{t-1} + \\epsilon_{CPI,t}$ and $IP_t = \\beta_{0,IP} + \\beta_{1,IP} IP_{t-1} + \\epsilon_{IP,t}$. In creating the lagged variables, you'll find the pandas command ```df.shift()``` useful.\n",
    "\n",
    "After you estimate ```mod = sm.ols(formula=\"CPI~CPI_1\", data=regr_df).fit()```, you can use ```mod.resid``` to access the error your forecast makes, that is \n",
    "\n",
    "$$F_{CPI,t} = CPI_t - (\\hat{\\beta}_{0,CPI} + \\hat{\\beta}_{1,CPI} CPI_{t-1})$$\n",
    "\n",
    "where the hat denotes the regression estimate.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### Problem 3b.\n",
    "Using your result from 3a and the excess returns of Apple, estimate the macroeconomic factor model. (Don't be surprised or disheartened by non-sensical or insignificant result. Macroeconomic factor models are known for this, and we didn't even build a sohpisticated one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Momentum\n",
    "\n",
    "The usual approach these days is to enhance the CAPM by the factor SMB (small minus big) and HML (high minus low). These factors are in your data file as well, and they are carefully constructed from firm characteristics.\n",
    "\n",
    "#### Problem 4a.\n",
    "\n",
    "Let's see what happens to our result for Apple and estimate\n",
    "\n",
    "$$R_{AAPL,t}^e = \\alpha_{AAPL} + \\beta_{AAPL, MKT} R_{MKT,t}^e + \\beta_{AAPL, SMB} SMB_t + \\beta_{AAPL, HML} HML_t +\\epsilon_{AAPL, t}.$$\n",
    "\n",
    "In particular, look at $R^2$ which measures how much variability of the stock returns the model explains, and the factor loading for the market factor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There has been some debate about a so-called momentum factor. Taking our idea from the autoregression of the macroeconomic factors, we want to look of the stock returns -- unexpected to our efficient market hypothesis -- have maybe a little bit of predictive content. The argument is that this tiny bit of content over time is very useful, although it does not offer any consistent effect in the short run.\n",
    "\n",
    "#### Problem 4b.\n",
    "\n",
    "Estimate \n",
    "\n",
    "$$R_{AAPL,t}^e = \\beta_{0,AAPL} + \\beta_{1,AAPL} R_{AAPL,t-1}^e + \\epsilon_{AAPL,t}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Outlook to the Homework.\n",
    "\n",
    "The momentum strategy builds on the observation that trends persist in some way or another. In your homework, you will therefore try to build a momentum factor (To make it a little more managable, we will restrict ourselves to the stocks of 10 companies, which you may choose at random. This will bias our results, but the procedure is what we're trying to understand.)\n",
    "\n",
    "If we get through the material faster, here are the steps outlined for you to try and go ahead:\n",
    "<ul>\n",
    "<li>Each month, rank the stocks by performance based on their return from best to worst.\n",
    "<li>Put the 30\\% best stock in a portfolio with weights $1/30$ each, and also include the 30\\% worst stock in there with weights $-1/30$ each. This means we go long the previous winners and sell short the previous losers.\n",
    "This portfolio is your factor! (In reality, we need to construct it more carefully, taking into account the presence of the factors SMB and HML, too.)\n",
    "<li>Run a regression for each stock including your new momentum factor, which will give you a bunch of $\\beta_{i,MOM}$.\n",
    "<li>Determine the average return of each asset, $\\langle R_{i,t}^e\\rangle.$ Now you have two lists: The average returns and the betas. Run a regression to find $\\lambda$.\n",
    "</ul>\n",
    "\n",
    "This procedure is called Fama-MacBeth regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
